Evolution Equations on Gabor Transforms and their Applications 



Remco Duits and Hartmut Fiihr and Bart Janssen and 
Mark Bruurmijn and Luc Florack and Hans van Assen 

October 28, 2011 

Abstract 

We introduce a systematic approach to the design, implementation and analysis of left-invariant 
evolution schemes acting on Gabor transform, primarily for applications in signal and image analysis. 
Within this approach we relate operators on signals to operators on Gabor transforms. In order to 
obtain a translation and modulation invariant operator on the space of signals, the corresponding 
operator on the reproducing kernel space of Gabor transforms must be left invariant, i.e. it should 
commute with the left regular action of the reduced Heisenberg group Hr- By using the left-invariant 
vector fields on Hr in the generators of our evolution equations on Gabor transforms, we naturally 
employ the essential group structure on the domain of a Gabor transform. Here we distinguish be- 
tween two tasks. Firstly, we consider non-linear adaptive left-invariant convection (reassignment) to 
sharpen Gabor transforms, while maintaining the original signal. Secondly, we consider signal en- 
hancement via left-invariant diffusion on the corresponding Gabor transform. We provide numerical 
experiments and analytical evidence for our methods and we consider an explicit medical imaging 
application. 

keywords: Evolution equations, Heisenberg group , Differential reassignment , Left-invariant vector 
fields , Diffusion on Lie groups , Gabor transforms , Medical imaging. 

1 Introduction 

The Gabor transform of a signal / G L2(M'^) is a function G^pif] : W'' x R'^ ^ C that can be roughly 
understood as a musical score of /, with G-,i>[f]{p, q) describing the contribution of frequency q to the 
behaviour of / near p [1, 2\. This interpretation is necessarily of limited precision, due to the various 
uncertainty principles, but it has nonetheless turned out to be a very rich source of mathematical theory 
as well as practical signal processing algorithms. 

The use of a window function for the Gabor transform results in a smooth, and to some extent blurred, 
time-frequency representation. For purposes of signal analysis, say for the extraction of instantaneous 
frequencies, various authors tried to improve the resolution of the Gabor transform, literally in order to 
sharpen the time-frequency picture of the signal; this type of procedure is often called "reassignment" in 
the literature. For instance, Kodera et al. |S] studied techniques for the enhancement of the spectrogram, 
i.e. the squared modulus of the short-time Fourier transform. Since the phase of the Gabor transform is 
neglected, the original signal is not easily recovered from the reassigned spectrogram. Since then, various 
authors developed reassignment methods that were intended to allow (approximate) signal recovery 
[HISJIS]. We claim that a proper treatment of phase may be understood as phase-covariance, rather than 
phase-invariance, as advocated previously. An illustration of this claim is contained in Figure [l] 

We adapt the group theoretical approach developed for the Euclidean motion groups in the recent 
works [SlliniiniinillllllllllllllSlin!, thus illustrating the scope of the methods devised for general 
Lie groups in (TS] in signal and image processing. Reassignment will be seen to be a special case of 
left-invariant convection. A useful source of ideas specific to Gabor analysis and reassignment was the 
paper [6^. 

1.1 Structure of the article 

This article provides a systematic approach to the design, implementation and analysis of left-invariant 
evolution schemes acting on Gabor transform, primarily for applications in signal and image analysis. 




Figure 1: Top row from left to right, (1) the Gabor transform of original signal /, (2) processed 
Gabor transform $t(yV^/) where $t denotes a phase invariant shift (for more elaborate adaptive con- 
vection/reassignment operators see Section |6] where we operationalize the theory in |6]) using a discrete 
Heisenberg group, where I represents discrete spatial shift and m denotes discrete local frequency, (3) 
processed Gabor transform $t(W^/) where <i>( denotes a phase covariant diffusion operator on Gabor 
transforms with stopping time t > 0. For details on phase covariant diffusions on Gabor transforms, 
see [21 ch:7] and "8^, ch:6]. Note that phase-covariance is preferable over phase invariance. For example, 
restoration of the old phase in the phase invariant shift creates noisy artificial patterns (middle image) 
in the phase of the transported strong responses in the Gabor domain. Bottom row, from left to right: 
(1) Original complex-valued signal /, (2) output signal T^/ = WJ^tW^/ where $t denotes a phase- 
invariant spatial shift (due to phase invariance the output signal looks bad and clearly phase invariant 
spatial shifts in the Gabor domain do not correspond to spatial shifts in the signal domain), (3) Output 
signal T^/ — W^'l'fW^/ where denotes phase- covariant adaptive diffusion in the Gabor domain with 
stopping time t > 0. 
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The article is structured as follows: 



Introduction and preliminaries: Sections [T] [2] In Section [2] we consider time-frequency anal- 
ysis and its inherit connection to the Heisenberg group and subsequently we introduce relevant 
structures on the Heisenberg group. 

• Evolution equations on the Heisenberg group and on related manifolds: Sections [3) 
|4]and[5j Section|3]we set up the left-invariant evolution equations on Gabor transforms. We explain 
the rationale behind these convection-diffusion schemes, and we comment on their interpretation 
in differential-geometric terms. Sections |4] and [5] are concerned with a transfer of the schemes from 
the full Heisenberg group to phase space, resulting in a dimension reduction that is beneficial for 
implementation. 

• Convection: In Section|6]we consider convection (reassignment) as an important special case. For 
a suitable choice of Gaussian window, it is possible to exploit Gauchy-Riemann equations for the 
analysis of the algorithms, and the design of more efficient alternatives. For example, in Section [6] 
we deduce from these Gauchy-Riemann relations that differential reassignment according to Daudet 
et al. [B], boils down to a convection system on Gabor transforms that is equivalent to erosion on 
the modulus, while preserving the phase. 

• Discretization and Implementation: Sections [7| and [Sj In order to derive various suitable 
algorithms for differential reassignment and diffusion we consider discrete Gabor transforms in 
Section [7j As these Gabor transforms are defined on discrete Heisenberg groups, we need left- 
invariant shift operators for the generators in our left-invariant evolutions on discrete Heisenberg 
groups. These discrete left-invariant shifts are derived in Section [8j They are crucial in the 
algorithms for left-invariant evolutions on Gabor transforms, as straightforward finite difference 
approximations of the continuous framework produce considerable errors due to phase oscillations 
in the Gabor domain. 

• Implementation and analysis of reassignment: Sections [9| and [TOl In Section |9] we employ 
the results from the previous Sections in four algorithms for discrete differential reassignment. We 
compare these four numerical algorithms by applying them to reassignment of a chirp signal. We 



provide evidence that it actually works as reassignment (via numerical experiments, subsection 9.1 ) 



and indeed yields a concentration about the expected curve in the time-frequency plane (Section 



10 1. We show this by deriving the corresponding analytic solutions of reassigned Gabor transforms 



of arbitrary chirp signals in Section [10) 

Diffusion: In Section [Tl] we consider signal enhancement via left-invariant diffusion on Gabor 
transforms. Here we do not apply thresholds on Gabor coefficients. Instead we use both spatial 
and frequency context of Gabor-atoms in locally adaptive flows in the Gabor domain. We include a 
basic experiment of enhancement of a noisy ID-chirp signal. This experiments is intended as a pre- 
liminary feasibility study to show possible benefits for various imaging- and signal applications. The 
benefits may be comparable to our directly relatecQ previous works on enhancement of (multiply 
crossing) elongated structures in 2D- and 3D medical images via nonlinear adaptive evolutions on 
invertible orientation scores and/or diffusion-weighted MRI images, [51 [TUl [TTl [HI [HI [Til [TSl [TBI [T7] . 

2D Imaging Applications: In Section [12) we investigate extensions of our algorithms to left- 
invariant evolutions on Gabor transforms of 2-dimensional greyscale images. We apply experiments 
of differential reassignment and texture enhancement on basic images. Finally, we consider a 
cardiac imaging application where cardiac wall deformations can be directly computed from robust 
frequency field estimations from Gabor transforms of MRI-tagging images. Our approach by Gabor 
transforms is inspired by [19 and it is relatively simple compared to existing approaches on cardiac 
deformation (or strain/velocity) estimations, cf. [201 [H] [HI [13) . 



^Replace H{2d-\-l) and the Schrodinger representation with SE{d) and its left-regular representation on L2(IR''), d = 2,3. 
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2 Gabor transforms and the reduced Heisenberg group 



Throughout the paper, we fix integers c? G N and n G Z \ {0}. The continuous Gabor-transform 
Gjf, [f] : M'' X M'' ^ C of a square integrable signal / : M'' ^ C is commonly defined as 



Gi. [/] (P, q)^ fim^- P)e-''^™ , (1) 



where G L2(M'') is a suitable window function. For window functions centered around zero both in 
space and frequency, the Gabor coefficient Q^f, [f] (p, q) expresses the contribution of the frequency nq to 
the behaviour of / near p. 

This interpretation is suggested by the Parseval formula associated to the Gabor transform, which 
reads 

/ / |e^[/](p,g)pdpd<? = C^ / |/(p)pdp, whereC^ = -|jV||^,(R.) (2) 
for all f,tp E L2(M'^). This property can be rephrased as an inversion formula: 

fiO = ^ j j ^v4/]fe9)e^'™'*"''^'V(^-p)dpdg , (3) 

R'i R"* 

to be read in the weak sense. The inversion formula is commonly understood as the decomposition of / 
into building blocks, indexed by a time and a frequency parameter; most applications of Gabor analysis 
are based on this heuristic interpretation. For many such applications, the phase of the Gabor transform 
is of secondary importance (see, e.g., the characterization of function spaces via Gabor coefficient decay 
[21] )■ However, since the Gabor transform uses highly oscillatory complex- valued functions, its phase 
information is often crucial, a fact that has been specifically acknowledged in the context of reassignment 
for Gabor transforms 

For this aspect of Gabor transform, as for many othert]^ the group-theoretic viewpoint becomes 
particularly beneficial. The underlying group is the reduced Heisenberg group H^.. As a set, H,. = 
M^'i X M/Z, with the group product 

(p, q, s + Z)(p', q', s'+Z) = {p + p', q + q',s + s' +^{q-p' -p- q') + Z) . 

This makes Hj. a connected (nonabelian) nilpotcnt Lie group. The Lie algebra is spanned by vectors 
Ai, . . . ,A2d+i with Lie brackets [^^,74^+^] = —A-id+i, and all other brackets vanishing. 
Hr acts on L2(M'^) via the Schrddinger representations U" : B{L.2{ 



The associated matrix coefficients are defined as 

q, s + Z) = {U^,,,,,^j^)i>, f)u{R-). (5) 

In the following, we will often omit the superscript n from U and W^, implicitly assuming that we use 
the same choice of n as in the definition of Q^. Then a simple comparison of ([5| with ([T]) reveals that 

e^[/](p,9) = W.0/(p,g,s = -^). (6) 

Since W^fip, q,s + Z) = e^''™''>V^/(p, g, + Z) , the phase variable s does not affect the modulus, and 
([2| can be rephrased as 

' |W^[/]fe9,5 + Z)pdpd<zds = C^ / |/(p)pdp. (7) 







^Regarding the phase factor that arises in the composition of time frequency shifts and the Hr group-structure in the 
Gabor domain we quote " This phase factor is absolutely essential for a deeper understanding of the mathematical structure 
of time frequency shifts, and it is the very reason for a non-commutative group in the analysis.'^ from 1241 Ch:9]. 
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Just as before, this induces a weak-sense inversion formula, which reads 

f I I W^[/]fe'?>5 + Z)Z^(p,5,.+z)^dpdgds . 

L^iA Jo JR<* JR'' 

As a byproduct of Q, we note that the Schrodinger representation is irreducible. Furthermore, the 
orthogonal projection of L2(-ffr) onto the range TZiW^) turns out to be right convolution with a 
suitable (reproducing) kernel function, 

{V^U){h) - U * K{h) = / U{g)K{g-^h)dg , 

with dg denoting the left Haar measure (which is just the Lebesgue measure on M^'* x E/Z) and 
K{p,q,s) = ^yV.4,tp{p,q,s) ^ ^(C/(p,9,s)'/'»- 

The chief reason for choosing the somewhat more redundant function W^/ over G-4>[f] is that 
translates time-frequency shifts acting on the signal / to shifts in the argument. If C and TZ denote the 
left and right regular representation, i.e., for all g,h E Hr and F e h2{Hr), 

{CgF){h) = F{g-'h) , {TlgF)(h) = F{hg) , 

then intertwines U and £, 

w^ou;; ^ Cgow-^ . (8) 

Thus the group parameter s in Hr keeps track of the phase shifts induced by the noncommutativity 
of time-frequency shifts. By contrast, right shifts on the Gabor transform corresponds to changing the 
window: 

T^g{W:^ih)) = [Uug^, /) = Wu^^Iih) . (9) 

3 Left Invariant Evolutions on Gabor Transforms 

We relate operators $ : TZ{W.^) — > L2(i/,) on Gabor transforms, which actually use and change the 
relevant phase information of a Gabor transform, in a well-posed manner to operators : L2(IR'*) — >■ 
L2(M'') on signals via 

(T^/)(O = (w;°<i'°w.0/)(C) 

= cr / / /($(Wv,/))(p,g,s)e''""[(«''')+"-(i/2^(P'«^V(C~p)dpdgds. (10) 
Our aim is to design operators that address signal processing problems such as denoising or detection. 
3.1 Design principles 

We now formulate a few desirable properties of T^, and sufficient conditions for $ to guarantee that T.^ 
meets these requirements. 

1. Covariance with respect to time- frequency- shifts: The operator T.0 should commute with time- 
frequency shifts; 

oUg ^Ug O T.jf, 

for all g E H{2d + 1). This requires a proper treatment of the phase. 

One easy way of guaranteeing covariance of T.^ is to ensure left invariance of <f>: 

<^> o Cg = £g o <i> (11) 

for all g G H{2d + 1). If $ commutes with Cg, for all g G H^, it follows from (jsj) that 

T^oui' = w; o $ o w^ou;; = w; o $ o £<, o w,p = u;' o t^, . 

Generally speaking, left invariance of 4> is not a necessary condition for invariance of T^: Note 
that Wl = WloV^. Thus if $ is left-invariant, and A : L2(Fr) -> TZiW^)-^ an arbitrary operator. 
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then $ + A cannot be expected to be left- invariant, but the resulting operator on the signal side 
will be the same as for $, thus covariant with respect to time-frequency shifts. 

The authors of [5] studied reassignment procedures that leave the phase invariant, whereas we 
shall put emphasis on phase- covariance. Note however that the two properties are not mutually 
exclusive; convection along equiphase lines fulfills both. (Sec also the discussion in subsection 3.4 ) 

2. Covariance with respect to rotation and translations : 

Tv,oZ^/^W=Z^^5^WoT^ (12) 

for all g G SE{d) = M'' xi SO{d) with unitary representation U^Eid) . g^f^^^ _^ B{h2{M.'^)) given 
^'^(x^R)' -/"fiO — /(^^^(C ~ 2;)), for almost every {x,R) € SE{d). Rigid body motions on signals 
and Gabor transforms relate via 

{Wu,,,^,^f){R-Hp-x),R-^q,s+if), 

for ah / e L2(M'^) and for all {x,R) G SE{d) and for all {p,q,s) e H{2d+ 1) and therefore we 
require the kernel to be isotropic (besides <& o £(2;, 0,0) — ^{x.0,0) ° ^ included in Eq. ( 11 )) and we 
require 

'j>°vj5? = v53fo<i> (14) 

for all ReSO{d). 

3. Nonlinearity: The requirement that T.^ commute with immediately rules out linear operators 

Recall that is irreducible, and by Schur's lemma [25], any linear intertwining operator is a 
scalar multiple of the identity operator. 

4. By contrast to left invariance, right invariance of $ is undesirable. By a similar argument as for 
left-invariance, it would provide that — Tj^n^. 

We stress that one cannot expect that the processed Gabor transform <i>(yV0/) is again the Gabor 
transform of a function constructed by the same kernel ip, i.e. we do not expect that $(7?.(>V^)) C 

3.2 Invariant differential operators on H^. 

The basic building blocks for the evolution equations are the left-invariant differential operators on Hr of 
degree one. These operators are conveniently obtained by differentiating the right regular representation, 
restricted to one-parameter subgroups through the generators {Ai, . . . , A2d+i \ = {dp^ , • • . , 9p^, , . . . , 9s} C 

d'R{A^)U{g) = lim ^^g^ ')-^(.9) ^ ^ smooth U £ C°°{Hr), (15) 

e-!-0 e 

The resulting differential operators {dTZ{Ai), . . . , dTZ{A2d+i)} ='■ {Ai, . . . , A2d+i} denote the left-invariant 
vector fields on Hr, and brief computation of (15 1 yields: 

At = -f ^9s, Ad+i = dq, - ^ds, A2d+i = ds, for i = 1, . . . , d. , 
The differential operators obey the same commutation relations as their Lie algebra counterparts Ai, . . . , A2d+i 

[Ai,Ad+i] ■— AiAd+i — Ad+iAi — ~A2d+i, (16) 
and all other commutators are zero. I.e. ATZ is a Lie algebra isomorphism. 
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3.3 Setting up the equations 

For the effective operator $, we will choose left-invariant evolution operators with stopping time t > 0. 
To stress the dependence on the stopping time we shall write <&( rather than $. Typically, such operators 
are defined by $t(yV^/)(p, q, s) = W{p, q, s, t) where W is the solution of 



dtW{p, q, S, t) = Q(|W^/|, ^1, . . . , A2d)W{p, q, s, t), 
W{p,q,s,0) ^W^f{p,q,s). 



(17) 



where we note that the left-invariant vector fields are given by 

Ai = dp^ + f 5s, Ad+i = dq^ - '-^ds,A2d+i ^ds, for i = 1, . . . , d, , 
with left-invariant quadratic differential form 

2d 2d 2d 

Q(| W^/l ,Ai,..., A2d) = -J2 «i(|W^/!)(P, l)Ai + Aid W^/|)(P, q) Aj. (18) 

i — 1 z — Ij — 1 

Here ai(|>V^/|) and Dij^W^f]) are functions such that {p,q) ai{\'yV^f\){p,q) G M and {p,q) i-> 
Dij{\yV^f\){p,q) G M are smooth and either D = (pure convection) or = D > holds pointwise 
(with D — [Dij]) for alH = 1, . . • , 2d, j = 1, . . . 2d. Moreover, in order to guarantee left-invariance, the 
mappings : W^/ i— >■ ai(|W^/|) need to fulfill the covariance relation 

a,i\ChW^f\){g) = a,i\W^f\){p ~p\q~ g'), (19) 

for all / e L2(M), and all g = {p, q, s + Z), /i = {p' , q' , s' + Z) e H^. 

For ai — ... — a2d+i — 0, the equation is a diffusion equation, whereas if = 0, the equation 
describes a convection. We note that existence, uniqueness and square-integrability of the solutions 
(and thus well-definedness of T) are issues that will have to be decided separately for each particular 
choice of a := {ai, . . . ,a2d)^ and D. In general existence and uniqueness are guaranteed, provided 
that the convection vector and the diffusion-matrix are smoothly depending on the initial condition, 
see Appendix [Xj Occasionally, we shall consider the case where the convection vector {ai)f^i and the 
diffusion-matrix D are updated with the absolute value of the current solution |Vl^(-,i)| at time t > 0, 
rather than having them prescribed by the modulus of the initial condition \W{-, 0)| = |W^(/)| = l^v/l 
at time t = 0, i.e. 

ai(|VF(-, 0)1) is sometimes replaced by ai(|VF(-, i)|) and/or 
Dij {\W{-,0)\) is sometimes replaced by Dij {\W{-,t)\) 

In case of such replacement the PDE gets non-linear and unique (weak) solutions are not a priori 
guaranteed. For example in the cases of differential re-assignment we shall consider in Chapter [6j we will 
restrict ourselves to viscosity solutions of the corresponding Hamilton- Jacobi evolution systems, |26[ I27j . 

In order to guarantee rotation covariance we set column vector a := (a^, . . . , a?'^)'^ and D = [Dij] 
with row-index i and column-index j and we require 

(a(vf„^ifl^(.,t)))(.9)=R(a(W^(,i)))(R-ig) , 
(i?(vj^^f P^(.,t)))(5) = R-^ (i?(m-,t)))(R-\g) R . 



for all R = e SO{2d), R G S'O(d), g G H{2d + 1), U G L2(i/(2d + 1)), where we recall ([131. 

This definition of $<, for each t > fixed, satisfies the criteria we set up above: 

1. Since the evolution equation is left-invariant (and provided uniqueness of the solutions), it follows 
that $t is left-invariant. Thus the associated is invariant under time-frequency shifts. 

2. The rotated left-invariant gradient transforms as follows 

UXo5?^(-'^)) =R {AR-igW{;t)){R-'g), withA={Al,...,A2df. 
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Thereby (the generator of) our diffusion operator $t is rotation covariant, i.e. 



Q{\W{-,t)\,A) o Vj^^f = VJ5J^ o Q{\W{-,t)\,A) for all R e SO{d), 



if Eq. (20 1 and Eq. il9n hold. For example, if a = and D would be constant, then by Eq. (20) and 



Schur's lemma one has D = diag{D", . . . , _D", . . . , Dd+i,d+i| yielding the Kohn's Laplacian 

= Dn Eti A' + Dd+i.d+i E'id+i -^l cf- EH], and indeed A;^ o vj^^f = Vfo^^i? o A^. 

3. In order to ensure non-linearity, not all of the functions a^, Dij should be constant, i.e. the schemes 
should be adaptive convection and/or adaptive diffusion, via adaptive choices of convection vectors 
(oi, . . . ,a2d)'^ and/or conductivity matrix D. We will use ideas similar to our previous work on 
adaptive diffusions on invertible orientation scores [22] , [H] , [H] , [H] . We use the absolute value of 
the (evolving) Gabor transform to adapt the diffusion and convection in order to avoid oscillations. 

4. The two-sided invariant differential operators of degree one correspond to the center of the Lie 
algebra, which is precisely the span of A2d+i- Both in the cases of diffusion and convection, we 
consistently removed the A2d+i — 9s-direction, and we removed the s-dependence in the coefficients 
ai{\y^i>f\){p,q), Dij{\W^f\){p,q) of the generator Q{\W^f\,Ai,...,A2d) by taking the absolute 
value |W^/|, which is independent of s. A more complete discussion of the role of the s- variable is 
contained in the following subsection. 

3.4 Convection and Diffusion along Horizontal Curves 



So far our motivation for (17) has been group theoretical. There is one issue we did not address yet. 



namely the omission of ds — A2d+i in (17 1. Here we first motivate this omission and then consider the 
differential geometrical consequence that (adaptive) convection and diffusion takes place along so-called 
horizontal curves. 

The reason for the removal of the A2d+i direction in our diffusions and convections is simply that this 
direction leads to a scalar multiplication operator mapping the space of Gabor transform to itself, since 
dsW^jjf = —2'KinW^f . Moreover, we adaptively steer the convections and diffusions by the modulus 
of a Gabor transform |W0/(p, g, s)| = \G)pf{p,q)\, which is independent of s, and clearly a vector field 
{p, q, s) !—> F{p, q)ds is left-invariant iff F is constant. Consequently it does not make sense to include 
the separate ds in our convection-diffusion equations, as it can only yield a scalar multiplication. Indeed, 
for all constant a > 0, /? € M we have 

[ds, Q{\W^f\,Ai, . . .,A2d)] = and dsW^^f = -2mnW^,f ^ 

gt((aO,"+/3a,)+Q(|W^/|,^i,...,^2d)) ^ ^~ta{27rnf-tl32^in gtQ(| W^/| ,^2d) _ 

In other words ds is a redundant direction in each tangent space Tg{Hr), g € Hr- This however does not 
imply that it is a redundant direction in the group manifold Hr itself, since clearly the s-axis represents 
the relevant phase and stores the non-commutative nature between position and frequency, [71 ch:l]. 

The omission of the redundant direction ds in T{Hj.) has an important geometrical consequence. 
Akin to our framework of linear evolutions on orientation scores, cf. [HI [22], this means that we enforce 
horizontal diffusion and convection. In other words, the generator of our evolutions will only include 
derivations within the horizontal part of the Lie algebra spanned by {Ai, A2}- On the Lie group Hr this 
means that transport and diffusion only takes place along so-called horizontal curves in Hr which are 
curves t {p{t),q{t), s{t)) e Hr, with s{t) e (0, 1), along which 

= l[ Y.q^{r)p[{r)-p^{r)q[{T)dT , (21) 

'=1 

see Theorem [l] This gives a nice geometric interpretation to the phase variable s{t), as by the Stokes 
theorem it represents the net surface area between a straight line connection between (p(0), ^(O), s(0)) and 
{p{t) , q(t) , s{t)) and the actual horizontal curve connection [0,t] 9 t i-> (p(t), g(r), s(r)). For example, 
horizontal diffusion with diagonal D is the forward Kolmogorov equation of Brownian motion t i— >■ 



{P{t), Q(t)) in position and frequency and Eq. (21 ) associates a random variable S{t) (measuring the net 
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surface area) to the implicit smoothing in the phase direction due to the commutator [^i, A2] = A3 = ds, 

cf. El HE]. 

In order to explain why the omission of the redundant direction ds from the tangent bundle T{Hr) 
implies a restriction to horizontal curves, we consider the dual frame associated to our frame of reference 
{Ai, . ■ . ,^2d+i}- We will denote this dual frame by {d^"'^, . . . ,dA?'^^^} and it is uniquely determined 
by {dA\Aj) = 5^ , i,i = 1,2,3 where (5j denotes the Kronecker delta. A brief computation yields 

AA'\ , ,= dp* , dA'^+'\ , , = dgS i = 1, . . . , d 

d-^''^'U,,,)-ds+i(p.dg-g.dp), ^''^ 
Consequently a smooth curve 1 1— > ^(t) = ^(t), •s(i)) is horizontal iff 

(d^^'^+ll^^^) ,7'(=s)) = ^ At) = hqit)-p\t)-pit) ■ q'{t)). 



Theorem 1 Let f G L2(M) be a signal and W^,f be its Gabor transform associated to the Schwartz 
function tp. If we just consider convection and no diffusion (i.e. D ^ 0) then the solution of is 
given by 

W{g, t) = W^/(7f (0) , .9 - <i, e ^r, 

where the characteristic horizontal curve t 1— >■ "fj"{t) ~ {p(t),q{t), s{t)) for each go — {po,qo,So) £ Hr is 
given by the unique solution of the following ODE: 

---a\\W.^f\)ip{t),q{t)), p(0)=po, 
---a^{\W^f\){p{t),q{t)), q(0) = go, 
= iffip(i) - mq^t), s(0) = so, 

Consequently, the operator W^f H' W{-,t) is phase covariant (the phase moves along with the charac- 
teristic curves of transport) : 

arg{W{g, t)} = arg{W^,f}i^j) for all t > 0. (23) 

Proof First we shaU show that 57^ = 7/(0 for .9 ^ Hr and aU t > and aU / e L2(M). To 

this end we note that both solutions are horizontal curves, i.e. 




[dA' 



where dA'^ — ds + 2(pdq — qdp). So it is sufficient to check whether the first two components of the curves 
coincide. Let g = {p, q, s) £ Hr and define pe : IR+ ^ K, 9e : ^ K and : M+ ^ M, : M+ ^ M by 

p,{t) := {dp,gru^_^f W> , Pgit) ■= (dp, 7^/ W) , 

qe{t) {dq,gYul_,f ' '^^^^^ ^ ^ / ^^^^ ' 

then it remains to be shown that pe + p = Pg and pe + q = qg. To this end we compute 

^^(i) = a\\W^,Ug-if\){p,it),q,it)) = a\\Cg-iW^f\){p,{t),q,{t)) = (I W.0/|) (pe W + 9e W + 9), 

so that we see that (pe + 9e + 9) satisfies the following ODE system: 

' f^ip + p,){t) ^a\\W^,f\){pe{t)+p,qe{t) + q), t > 0, 
|(g + ge)W =a2(|W^/|)(ge(i) + 9,9eW+9), t>0, 

P+Pe=P, 

q + qe^q. 

This initial value problem has a unique smooth solution, so indeed Pg = p + Pe and qg = q + q^. Finally, 
we have by means of the chain-rule for differentiation: 

iiy^^f){l%t)) =E(dAl,.(t),7fW> (A|,.(,)W^/)(7K0) 

= M^^t)^^f{i'f(t)) + i9(t) ^2|^.(,)W^/(7?(i)) 

= -ai(|W>>/l)bsW,95W)^i|^«(t)W^(7fW) 

-a\\W^f\){pg{t),qg{t)) ^2|^.(,) W^(7fW) • 
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from which the result follows □ 

Also for the (degenerate) diffusion case with D = = [Dyli j=i,...,ti > 0, the omission of the 2d+l^^ 
direction dg — A2d+i implies that diffusion takes place along horizontal curves. Moreover, the omission 



does not affect the smoothness and uniqueness of the solutions of (17), since the initial condition is 
infinitely differentiable (if t/i is a Schwarz function) and the Horniander condition [3D], [TH] is by (16) still 
satisfied. 

The removal of the ds direction from the tangent space does not imply that one can entirely ignore the 
s-axis in the domain of a (processed) Gabor transform. The domain of a (processed) Gabor transform 
^•((W^/) should no^he considered as K^'^ = H^/Q. Simply, because [dp, dq] = whereas we should have 
( jlG] ). For further differential geometrical details see the appendices of [71, analogous to the differential 
geometry on orientation scores, [H], [3 App. D , App. C.l ]. 



4 Towards Phase Space and Back 

As pointed out in the introduction it is very important to keep track of the phase variable s > 0. The 
first concern that arises here is whether this results in slower algorithms. In this section we will show 
that this is not the case. As we will explain next, one can use an invertible mapping S from the space 
T-Ln of Gabor transforms to phase space (the space of Gabor transforms restricted to the plane s — 
As a result by means of conjugation with S we can map our diffusions on C L2(K^ x [0, 1]) uniquely 
to diffusions on L2(M^) simply by conjugation with S. From a geometrical point of view it is easier 
to consider the diffusions on T-Ln C L2(M^'^ x [0,1]) than on L2(M^'^), even though all our numerical 
PDF- Algorithms take place in phase space in order to gain speed. 

Definition 2 Let Hn denote the space of all complex-valued functions F on Hr such that F{p, q, s+Z) ~ 
e-^^'^->Fip,q,0) and F(-,-,s + Z) G h2{R'^'^) for all s € M. Clearly W^f G for all f,ij€nn. 

In fact Hn is the closure of the space {W^V \ ip, f £ L2(M)} in L2(-ffr)- The space Hn is bi-invariant. 
Since 

o = C, o w; and W^„^ = 7^, o , (24) 

where TZ denotes the right regular representation on L2(i/, ) and £ denotes the left regular representation 
of Hr on h2{Hr). We can identify Tin with L2(M^'^) by means of the following operator S : Hn — >■ L2(]R^'^) 
given by 

iSF)ip, q) - F{p, q,^+Z) = e™F(p, q,0 + Z). (25) 
Clearly, this operator is invertible and its inverse is given by 

iS-^F){p, q,S + Z)^ e'2^^sn^-^nnpqp^p^ ^26) 

The operator S simply corresponds to taking the section s{p,q) = — ^ in the left cosets Hr/Q where 
= {(0, 0, s + Z) I s G K} of Hr- Furthermore we recall the common Gabor transform Q"^ given by (jT| 
and its relation (joj) to the full Gabor transform, which we can now write as — S o W^. 

Theorem 3 Let the operator $ map the closure T-Ln, n g Z, of the space of Gabor transforms into itself, 
i.e. $ : "Hn — > 'Hn- Define the left and right-regular rep's of Hr on Hn by restriction 

T^gln,, <^rid = Cg\^^^ for all g e Hr. (27) 

Define the corresponding left and right-regular rep 's of Hr on phase space by 

n^^^ := s o 7^("^ o s-\ 4"^ := s o o s-\ 

For explicit formulas see /?, p. 9]. Let $ 5 o $ o S^^ be the corresponding operator on L2(1R^'') and 



^As we explain in [3 App. B and App. C ] the Gabor domain is a principal fiber bundle Pj- = (H^, T, tt, 7?.) equipped 
with the Cartan connection form Ug (Xg ) = (ds + ^ (pdg — qdp) , Xg), or equivalently, it is a contact manifold, cf. [311 p. 6] , 
[3 App. B, def. B.14] , (Hr,d^2d+i)) 
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Then one has the following correspondence: 



(28) 



If moreover <^{TZ{yV^)) C TZiW^p) then the left implication may he replaced by an equivalence. If ^ does 
not satisfy this property then one may replace $ — >■ W^W^^ in (28) to obtain full equivalence. Note 

that = w;$w^ = w;(w^w;$)w^. 

Proof For details see our technical report [?! Thm 2.2]. 



5 Left-invariant Evolutions on Phase Space 



For the next three chapters, for the sake of simplicity, we fix d = 1 and consider left-invariant evolutions 
on Gabor transforms of ID-signals. We will return to the case d = 2 in Section [12] where besides of some 
extra bookkeeping the rotation-covariance {2nd design principle) comes into play. 

We want to apply Theorem|3]to our left-invariant evolutions ( 17) to obtain the left-invariant diffusions 
on phase space (where we reduce 1 dimension in the domain). To this end we first compute the left- 
invariant vector fields {.4^} := {SAiS^^}^^^ on phase space. The left-invariant vector fields on phase 
space are 

AiU{p' ,q') = SAxS-^U{p' ,q') = {{dp^~2n-Kiq')U){p',q'), 

A2U{p',q')=SA2S-'U{p',q') = {d,,U){p',q'), (29) 
A3U(p',q')=SA3S-'U{p',q') = -2innU{p',q') , 

for all (p, (7) e M and all locally defined smooth functions U : ^{p^q) C ^ C. 

Now that we have computed the left-invariant vector fields on phase space, we can express our left- 



invariant evolution equations (17) on phase space 



dtW{p, q, t) = Q{\g.^f\,Ai,A2)W{p, q, t), 
Wip,q,0) =g^.f{p,q). 



with left-invariant quadratic differential form 

Q{\g^f\,Ai,A2) 



E««(ie*/l)(p,9)A + ;^^A D,,{\g^f\){p,q)A, 

1—1 i—l j—1 



(30) 



(31) 



Similar to the group case, the and Dij are functions such that (j),q) 1— t- ai{\Q^f\){p,q) e M and 
{p,q) 1-^ o,i{\Q.^f\){p,q) e M are smooth and either D = (pure convection) or — D > (with 
D — [Dij] i,j = l,...,2d), so Hormander's condition |30j (which guarantees smooth solutions W, 
provided the initial condition W{-,-,0) is smooth) is satisfied because of (16). 



Theorem 4 The unique solution W of (30) is obtained from the unique solution W of (11) by means of 

W{p,q,t) = {SW{-rr,t)){p,q) , for all t > and for all {p,q) G M^ 
with in particular W{p,q,0) ^Q^{p,q) = {S'W^){p,q) = {SW{-,-,-,Q)){p,q) . 



Proof This follows by the fact that the evolutions (17) leave the function space "Hn invariant and the 
fact that the evolutions (30) leave the space L2(M^) invariant, so that we can apply direct conjugation 



with the invertible operator S to relate the unique solutions, where we have 
W{p, q, t) = (e*'3(ie^/l^^i-^=)g^/)(p, q) 

^ (g5otQ(|Wv,/|,-4i,^2)o5-i5y|;^j)(p^^) 

= (5 o e'Q(mJl^i-A2) o S-^SoW4,f){p,q) 

^{SW{;;;mp,q) 



(32) 



for all t > on densely defined domains. For every ?/; e L2(M) n 5'(M), the space of Gabor transforms is 
a reproducing kernel space with a bounded and smooth reproducing kernel, so that W^/ (and thereby 
|yV^/| = 1^/^/1) is uniformly bounded and continuous and equality (32) holds for all p,q G M^. □ 
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5.1 The Cauchy Riemann Equations on Gabor Transforms and the Under- 
lying Differential Geometry 

As previously observed in [6], the Gabor transforms associated to Gaussian windows obey Cauchy- 
Riemann equations which are particularly useful for the analysis of convection schemes, as well as for 
the design of more efficient algorithms. More precisely, if the window is a Gaussian tp{^) = i^aiO ■= 

e and / is some arbitrary signal in L2(M) then we have 



(a-Ua + iaAi)W^if) = ^ (a^Ua + iaAi) logW^if) = (33) 

where we included a general scaling a > 0. On phase space this boils down to 

{a-^A2+iaAi)g4,{f) = (34) 

since g^{f) = 5W^(/) and A, = S-^A,S for i = 1,2, 3. 

For the case a = 1, equation (34) was noted in [B]. Regarding general scale a > we note that 

G->Paf{P,l) = y/a gvi4,{f){p,q) = y/ag^,Vif{^,aq) 
with ■(/; = ipa=i with unitary dilation operator 2?q : L2(M) — > L2(M) given by 

Va{i^){x)=a-if{x/a), a>0. (35) 



As a direct consequence of Eq. ( 34 1 , respectively ( 33 ) , we have 

iLf^ldqCl'' ^ -a^dplU^l and \U''\dpn'' = a-^dqlU^l + 2TTq. 

A,n'^ = -a^^ and A^n'^^a'^^. ^^^> 

where C7" = G^i-Af), U^= WVa(/), ^^" = arg{e^„(/)} and l^'' = arg{W^„ (/)}. 

The proper differential-geometric context for the analysis of the evolution equations (and in particular 
the involved Cauchy-Riemann equations) that we study in this paper is provided by sub-Riemannian ge- 
ometry. As already mentioned in Subsection 3.4 we omit A3 from the tangent bundle T[Hj.) and consider 
the sub-Riemannian manifold {Hr,dA'^), recall (22), as the domain of the evolved Gabor transforms. 
Akin to our previous work on parabolic evolutions on orientation scores (defined on the sub-Riemannian 
manifold {SE{2), — sin0da; + cos6'dy)) cf.[T^, we need a left-invariant first fundamental form (i.e. metric 
tensor) on this sub-Riemannian manifold in order to analyze our parabolic evolutions from a geometric 
viewpoint. 

Lemma 5 The only left-invariant metric tensors : Hr x T{Hr) x T{Hr) C on the sub-Riemannian 
manifold {Hr,dA'^) are given by 

{i.j)€{l,2}^ 

Proof Let : i/,. x T{Hr) x T{Hr) — C be a left-invariant metric tensor on the sub-Riemannian 
manifold {H,.,dA^). Then since the tangent space of {Hr,dA'^) at g G Hr is spanned by {-^il^ , y^al^} 
we have 

®9= 9^J{9)dA'\^(g> dA^g 

for some gij{g) G C. Now & is left-invariant, meaning &gfi{{Lg)^Xh, {Lg)^Yh) — <3hiXii,Yii) for all 
vector fields X,Y on Hr, and since our basis of left-invariant vector fields satisfies (ig)* Ai\i^ = -^ilgh 
and {dA^,Aj) = (5* we deduce gijigh) — gij{h) = gij{e) for all g, h £ Hr- □ 

Akin to the related quadratic forms in the generator of our parabolic evolutions we restrict ourselves to 
the case where the metric tensor is diagonal 

&p= 9ij'^^' ® dy^^ = /^^dy^i (g) dA^ + dA'^ ® dA"^. (37) 

(»J)6{1,2}2 
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Here the fundamental positive parameter /3~^ has physical dimension length, so that this first funda- 
mental form is consistent with respect to physical dimensions. Intuitively, the parameter /3 sets a global 
balance between changes in frequency space and changes in position space within the induced metric 

7(0) = 9,7(1) = h, 

Note that the metric tensor &j3 is bijectively related to the linear operator G : — > 55', where S) ~ 

spanjy^i, ^^2} denotes the horizontal part of the tangent space, that maps Ai to and A2 to d.4^. 
The inverse operator of G is bijectively related to 

0^^ = ^ g'^At «) Aj = l3~'^Ai (»Ai+A2<»A2 ■ 

The fundamental parameter /3 is inevitable when dealing with flows in the Gabor domain, since position 
p and frequency q have different physical dimension. Later, in Section we will need this metric in 
the design of adaptive diffusion on Gabor transforms. In this section we primarily use it for geometric 
understanding of the Cauchy-Riemann relations on Gabor transforms, which we will employ in our 
convection schemes in the subsequent section. 

In potential theory and fluid dynamics the Cauchy-Riemann equations for complex-valued analytic 
functions, impose orthogonality between flowlines and equipotential lines. A similar geometrical inter- 



pretation can be deduced from the Cauchy-Riemann relations ( 36 1 on Gabor transforms : 
Lemma 6 Let U := W^^/ = |C/|e*^ be the Gabor transform of a signal f E L2(M). Then 

©-^^(dlog|f/|,Py,.df]) = ©-^i(d|t/|,P.^-df^) = 0, (38) 

where the left-invariant gradient of modulus and phase equal 

3 2 
dO = ^Af^dX\ d|C/| =^A|C/| dy^\ 

i=l i=l 
2 

whose horizontal part equals Vs^*dfl — ^ Aiil dA^ , Vf,td\U\ = d\U\. 



Proof By the second line in ( 36 1 we have 



« (5^^„-i(dlog|i7|,Pi5.dl7) = a — +a — = 0, 

from which the result follows. □ 

Corollary 7 Let go G Hr- Let W^,J = |C/|e'" with in particular Wv-„/(5o) = |C/(5o)|e*"'^"-'- Then the 
horizontal part Psj" 'i-^lgg of the normal covector dJlj^^ to the equi-phase surface {{p, q, s) € Hr \ ^{p, q, s) — 
f2(go)} is <&p- orthogonal to the normal covector d\U\\g to the equi- amplitude surface {(p, s) £ \ \U\{p^ q, s) 

\U\igo)}- 

For a visualization of the Cauchy-Riemannian geometry see Fig[2] where we also include the exponential 



curves along which our diffusion and convection (Eq.(17l) take place 



Remark 8 Akin to our framework of left-invariant evolutions on orientation scores fT^ one express 
the left-invariant evolutions on Gabor transforms (Eq. ^11^ ) in covariant derivatives, so that transport 
and diffusion takes place along the covariantly constant curves (auto-parallels) w.r.t. Cartan Connection 
on the sub-Riemannian manifold [Hr,dA^). A brief computation (for analogous details see shows 
that the auto-parallels t i-^ ^{t) w.r.t. the Cartan connection V coincide with the horizontal exponential 
curves. Auto-parallels are by definition curves that satisfy 

y^i = 0^f-J2 4,7'7^' = f = 0, (39) 
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where in case of the Cartan- connection the Christoff el- symbols coincide wit h m inus the anti- symmetric 
Lie algebra structure constants and with 7' — (d,4*| ,7). So indeed Eq. (39) holds iff Y = c* S K, 

2 

i = 1, 2, i.e. 7(t) = 7(0) exp{t J2 c'^i) for all t e R. 




Figure 2: Equi-aniplitude plane (red) and equi-phase plane in the Gabor transform of a chirp signal, 



that we shall compute exactly in Section [TO 



The left-invariant horizontal gradients of these surfaces are 
~ The 



locally ©^-orthogonal to each other, cf. Eq. (38), with 0^ = f^'^dA^(^dA'^+dA'^(^dA'^ with (3 = 



left-invariant vector fields {^1, y^2, -^3} = {dp + |(3s, dq — |9s, ds} form a local frame of reference which 
is indicated by the arrows. Some exponential curves (the auto-parallel curves w.r.t. Cartan connection) 
are indicated by dashed lines. 



6 Convection operators on Gabor Transforms that are both 
phase-covariant and phase-invariant 

In differential reassignment, cf. [SI [5] the practical goal is to sharpen Gabor distributions towards lines 
(minimal energy curves [71 App.D]) in Hr, while maintaining the signal as much as possible. 

We would like to achieve this by left-invariant convection on Gabor transforms U := yV^(/). This 



means one should set D = Oin Eq. (171 and ( 30 1 while considering the mapping U t-^ W{-,t) for a suitably 
chosen fixed time t > 0. Let us denote this mapping by $t : Hn — > H„ given by ^t{U) = W{-,t). Such 
a mapping is called phase invariant if 

arg (($t(t/))(p, q, s)) = arg (C/(p, q, s)) , 

for ah {p,q,s) e Hr and ah U G Hn, allowing us to write $t(|t7|e*") = e*" • <i>f''WU\) where is the 
effective operator on the modulus. Such a mapping is called phase covariant if the phase moves along 



with the flow (characteristic curves of transport), i.e. if Eq. (23) is satisfied. Somewhat contrary to 
intuition, the two properties are not exclusive. 



Our convection operators (obtained by setting D = in Eq. (17)) are both phase covariant and 
phase invariant iff their generator is. In order to achieve both phase invariance and phase covariance 
one should construct the generator such that the flow is along equi-phase planes of the initial Gabor 
transform W^f. As we restricted ourselves to horizontal convection there is only one direction in the 
horizontal part of the tangent space we can use. 
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Lemma 9 let rig^ be an open set in Hr and let U : ^ 



90 



be dijferentiable. The only horizontal 



direction in the tangent bundle above flg^ that preserves the phase of U is given by 

span{-A2^Ai + AinA2}, 

with = arg{t/}. 

Proof The horizontal part of the tangent space is spanned by {y^i, .42}, the horizontal part of the phase 
gradient is given by Pj^.dr? = y^i^ldy^-^ + A2^dA^. Solving for 



□ 



{F^ydn,a^Ai +a^A2) = 

yields = A(-^2f^, -4if)), A e C. 

As a result it is natural to consider the following class of convection operators. 

Lemma 10 The horizontal, left- invariant, convection generators C : ~^ T~in given by 

C{U) = M{\U\){-A2^AiU + AinA2U), where Q = arg{f7}, 

and where A4{\U\) a multiplication operator naturally associated to a bounded monotonically increasing 
differentiable function fi : [0, max(|;7|)] [0, ^(max(|C/|))] C M with ^{0) = 0, i.e. (A^(|C/|)F)(p, g) = 
/i(|J7|(p, (7)) V{j),ci) for all V G 'H„,(p, G M^, are well-defined and both phase covariant and phase 
invariant. 

Proof Phase covariance and phase invariance follows by Lemma |9] The operators are well-posed as the 
absolute value of Gabor transform is almost everywhere smooth (if ^ is a Schwarz function) bounded 
and moreover C can be considered as an unbounded operator from Hn into Hn, as the bi-invariant space 
"Hn is invariant under bounded multiplication operators that do not depend on the phase z = e"^^^^. □ 



For Gaussian kernels ipaiO — ^ ^ may apply the Cauchy Riemann relations (34 1 which simpli- 

fies for the special case A^(|C/|) — \U\ to 



C{e'''\U\) = {a'{dp\U\y + a-'id,\U\y) e^'K 
Now consider the following phase-invariant adaptive convection equation on H^. 



dtWig,t)^-CiWi;t)){g), 
Wig,0) = Uig) 



with either 



1. CiW{;t)) ^ M{\U\) i-A2n,Aifl) ■ {AiW{;t),A2Wi-,t)) or 



2. C{W{;t)) 



■in {„2{d^\W{-,t)\f 



\W(;t)\ 



-2 (d,\W{-,t)\)' 
\W{-.t)\ 



(40) 



(41) 



(42) 



In the first choice we stress that arg(M^(-,t)) = arg(M^(-,0)) = fi, since transport only takes place along 
iso-phase surfaces. Initially, in case A^(|C/|) = 1 the two approaches are the same since at i = the 
Gauchy Riemann relations ( 36 ) hold, but as time increases the Gauchy-Riemann equations are violated 



(this directly follows by the preservation of phase and non-preservation of amplitude). Gonsequently, 
generalizing the single step convection schemes in [HI [5] to a continuous time axis produces two options: 



1. With respect to the first choice in (42) in (|4ip (which is much more cumbersome to implement) we 



follow the authors in [B] and consider the equivalent equation on phase space: 
dtW{p,q,t)^~CiW{;t)){p,q), 

W{p,q,0)^g^f{p,q) -: C7(p, q) = e»f^(f |(7(p, g) | = e'f^(f |(7| (p, g) 



(43) 



with C{W{-,t)) = M{\U\) {^-A2nAiWi-,t) + idgn-2TTq)A2W{-,t)j, where we recah = SW^ 

and Ai = S^^AiS for i = 1,2, 3. Note that the authors in [6 consider the case A4 = 1. In addition 
to [5] we provide in Section |9] an explicit computational finite difference scheme acting on discrete 
subgroups of Hr, which is non-trivial due to the oscillations in the Gabor domain. 



15 



The second choice m Eq. (42) within Eq. (41 1 is just a phase-invariant inverse Hamilton Jacobi 
equation on H^, with a Gabor transform as initial solution. Rather than computing the viscosity 
solution cf. |26j of this non-linear PDE, we may as well store the phase and apply an inverse 
Hamilton Jacobi system on with the amplitude \U\ as initial condition and multiply with the 
stored phase factor afterwards. More precisely, the viscosity solution of the Hamilton Jacobi system 
on the modulus is given by a basic inverse convolution over the (max, -|-) algebra, [32) . (also known 
as erosion operator in image analysis) 



W{p, q, t) = ($t({/))(p, q, t) = [Kt e |C/|)(p, g)e'"(f ^'■*) , (44) 



KAV.,)^"^^^^^ (45) 



(/eg)(p,g)= inf Jff(p',(z') + /(P-P','z-'?')]- 



with kernel 



where 



Here the homomorphism between erosion and inverse diffusion is given by the Cramer transform 
C = 5^ologo£, [22], [33] I that is a concatenation of the multi-variate Laplace transform, logarithm 
and Fenchel transform so that 

c(/ * <?) = 5iog/:(/ * 5) = ^(log/:/ + log £5) = c/ e Cg, 

with convolution on the (max, -|-)-algebra given by / ® 3(x) = sup [/(x — y) + .9(y)]- 



7 Discrete Gabor Transforms 

In order to derive various suitable algorithms for differential reassignment and diffusion we consider 
discrete Gabor transforms. We show that Gabor transforms are defined on a (finite) group quotient 
within the discrete Heisenberg group. In the subsequent section we shall construct left-invariant shift 
operators for the generators in our left-invariant evolutions on this quotient group. 

Let the discrete signal is given by f = {f["]}^=o^ := {/ {-^)}n=a ^ ■ Let the discrete kernel be 
given by a sampled Gaussian kernel 

(|„|_L«zlij)2^ 

with ^ — where is the variance of the Gaussian. The discrete Gabor transform of f is then given 

by 

{WRi)[l,m,k] e-^''*^*-^)^ J2 xP[n - lL]i[n]e-^^ (47) 
with L, K, N, M,Q eN and 

N 

fc = 0, 1,...,Q-1 and/ = 0, ...,^^-1, m = 0, ...,M-1, with L = — , (48) 

K 

and integer oversampling P = M/L E Z. Note that we follow the notational conventions of the review 
paper |34 . One has 

- = (49) 
M PN- ^ ' 

It is important that the discrete kernel is N periodic since N = KL implies 

yfe^M)^Ln^,kW^nl + K,m,k]=W^^[l,m,k]^y^=o.,,,N^p[n-N]=^P[n], 



where / — {0, . . . , iV — 1}. Moreover, we note that the kernel chosen in (46) is even. 

For Riemann-integrable / with support within [0, 1] and "0 even with support within [—1, 1], say 

'/'(0 = e-^^^l[_ia](0, (50) 
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we have 





Consequently, we obtain the pointwise hmit (in the reproducing kernel space of Gabor transforms) 

/ tn f\ ft* 

(W-^f)[/, m, A:] ^ = -, g = — , s = -) as TV -> oo, (51) 

where we keep both P and -ftT fixed so that only M — oo as — ?► oo and with with scaled Gaussian 
kernel V(0 = e"''°~ * To this end we recall that the continuous Gabor transform was given 

by 

WrV(p,9,5)-e-2"(-¥) / ^(C-p)/(e)e-2"«'dC. 



7.1 Diagonalization of the Gabor transform 

In our algorithms, we follow [35j and [34j and use the diagonalization of the discrete Gabor transform 
by means of the discrete Zak-transform. The finite frame operator ^ : i2{I) ^2(1) equals 

K-l A/-1 

[^flM = E E('^'™'f)^'™N, ne/ = {0,...,iV-l}, 

1=0 m=0 

with xpi^ = „j Qim ^xl). Operator = i?* is coercive and has the following orthonormal eigenvectors: 

00 

Unk[n'] = -^v[n' ~ n] ("'""), with v{n) = V <5[n - IL] 

/— — 00 

for n e {0, . . . , L — 1}, fc e {0, . . . , if — 1}, where N = KL, and it is diagonahzed by: 

L-l A'-l 

with Discrete Zak transform given by [Z^f|[n,fc] — {Unk,i)i2{i)j i-^- — -^ri/£(urifc, f)u„fc, with 

p-1 

eigenvalues A„fc — L \Zxp^[n, k — PjjW^ and integer oversampling factor P = M/L. 

8 Discrete Left-invariant vector fields 

Let the discrete signal is given by f = {fMI^^Jo^ '■— {/ {-^)}n=Q ^ ■ Then similar to the continuous 
case the discrete Gabor transform of f can be written 

[>Vjfl[;,m,fc] = (^,,„,fc]t/>,f),,(,), (52) 

JV-l 

where / = {0, . . . , iV — 1} and (a, b) = ^ ^ Uibi and where 

Uii^mM^n] = e2"(t-S^)e'^^[„ _ n], (53) 
I — 0, . . . , K — l,m = 0, . . . , M — 1, k = 0, . . . ,Q — 1. Next we will show that (under minor additional 



conditions) Eq. (53) gives rise to a group representation of a finite dimensional Heisenberg group t)r, 



obtained by taking the quotient of the discrete Heisenberg group hr with a normal subgroup. 
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Definition 11 Assume ^ € N then the group hr is the set Z'^/({0} x {0} x QZ) endowed with group 
product 



[l,m,k][l',m',k'] = [l + l',m + m',k + k' + ^{ml' - m'l)ModQ]. 



(54) 



Lemma 12 Assume ^ G N and ^ G N and L even, N even. Then the subgroup 

[KZ, MZ, QZ] := {[IK, mM, kQ] \l,m,keZ} 
is a normal subgroup of hr. Thereby, the quotient t)r hr/{[KZ, MZ,QZ]) is a group with product 

[l,m,k][l',m',k'] = [l + l' ModK,m + m' ModM,k + k' + ^{ml' -ml) Mod{Q)]. (55) 



Eq. (53) defines a group representation on this group and thereby \)r is the domain of discrete Gabor 
transforms endowed with the group product given by Eq. [55). 

Proof Direct computation yields 

[I, m, k][l'K, m'M, k'Q] = [I + I'K, m + m'M, k + k'Q + ^{ml'K - m! Ml) + ModQ] = 
\!K, m'M, k'Q] [I, m, k] = [I + I'K, m + m'M, k + k'Q~ ^ {ml'K - m'Ml) + ModQ] 

since M/P = L G N and we assumed K/P G N. Consequently, H := [KZ,MZ,QZ] is a normal 
subgroup and thereby (since giHg2H] ~ gig2H) the quotient f)r ■— hr / {[KZ, MZ, QZ]) is a group with 
well-defined group product (55 1. The remainder now follows by direct verification of 



^[l,m,k]^ll',m',k'] — U[i^rnM][l' ,rn' ,k'] 



and „j fc] —U^l + K,m+M,k+Q]- 



□ 



In view of Eq. (51) we define a monomorphism between hr and Hr as follows. 



Lemma 13 Define the mapping (p : hr ^ Hr by 



4>[l, m, k] 



I mK k 



K' P ' Q, 

which sets a monomorphism between hrand the continuous Heisenberg group Hr 
Proof Straightforward computation yields 



(j)[l,m,k]<j)[l',m',k'] = 



J_ mK 
K ' P 



4)( 



/ mK k_ 
K' P ^ Q 



(^IM^^ irn+rn')K ^ k+k' + ^ml' -Im') ^^ = m, fc] [/', to', fc'] 



from which the result follows. □ 
The mapping (/) maps the discrete variables on a uniform grid in the continuous domain: 

s G [0,1) fc G [0,Q) nz, pe[0,l) ^ I e[0,K)nZ, g G [0, iV) m G [0, Af ) n Z. 

On the quotient-group t)r we defi ne th e for ward left-invariant vector fields on discrete Gabor-transforms 
as follows (where we again use ( 49 1 and ( 48 ) ) : 



(A?^ W^m, m,k] (d7^°+ [1, 0, 0]W^m, m, k] 



e P W^i[l + l,m,k]-Wj^t[l,m.k] 



{Atw^m,m,k] = M(d7^^+[o,l,o]Vy^f)[/,m,fc] = 



[Afw^r)[l,m,k] 



= Q(d7^^ [0,0,l]W^^f)[Z,m,fc] = 



W^f([l,m,k\[lfi,0\)-W^t[l,m,k] 



e M wj^{[l + 1.7n,k]-W^{[l,m,k] 

e+T^ VK^f[i,m + l,fc]-W^f[!,m,fc] 
K p-1 ' 

e+TT W^f[I,m + l,fe]-VK^f[I,m,fcl 
JV ' 

Wj^t[l.m,k + l]~wj^t[l,m.k] 

Q[e=^ -l)W^f[l,m,k]), 



(56) 
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and the backward discrete left-invariant vector fields 

VK'? f[/,m,fc]-e+T^ W^f[(-l,m,/c] 

{A?-W{^m,m,k] = {dn^-[l,0,0]W^m,m,k] = j^^^ , 

{ArW^m,^,k] = (d7^^+[0,l,0]Ty^f)[^,m,fc] = , 

iArW^m,m,k] = (d7^^+[0,0,l]W^^f)[^,m,fc] = ^ =Q{l-e^)W^f[l,m,k] 

(57) 



Remark 14 With respect to the step-sizes in (56) and (51) we have setp — j^, Q — "57-^7 ^ ~ W ' ^ ~ q ' 
so that the actual discrete steps are Ap — , Aq = N and As = Q^^ . This discretization is 
chosen such that both the continuous Gabor transform and the continuous left- invariant vector fields 



follow from their discrete counterparts by N 00, e.g. recall Eq. (51) 



Akin to the continuous case we use the following discrete version S of the operator S that maps a 
Gabor transform onto its phase space representation G^: 

The inverse is given by W^i[l,m, k] = {{S'^)-^G'^f)[l,m,k] = e"^"(^ + w)^^ fj^^ 

Again we can use the conjugation with to map the left-invariant discrete vector fields 
to the corresponding discrete vector fields on the discrete phase space: Af = {S^) o Af o (5^)^^. A 
brief computation yields the following forward left-invariant differences 

{ArG^m.m] = 

{AfG^m^m] ^MN-UG^f[l,m + l]~G^{[l,m]) (58) 
{AfG^t)[Lm] ^Q{e=^ -l)G^{[l,m] 

and the following backward left-invariant differences: 

(Gf, f)[/,m]-eT^G° f[!-l,m] 

(^f-G^f)[?,ml ^ 



[A^-G^mM =A/iV-i(Gjf[Z,m]-Gjf[/,™-l]) (59) 
{A^-G^mM ^Q{l-e'^)G^i[lM- 

The discrete operators are defined on the discrete quotient group fi^ and do not involve approximations 
in the setting of discrete Gabor transforms. They are first order approximation of the corresponding 
continuous operators on Hr as we motivate next. For / compactly supported on [0, 1] and both / and '0 
Riemann-integrable on M: 



-e^/^'(^-^)/(e)e-^^de 

K 

A^— 1 

O(^) - ^ V^' (li - ^) / (IL) e-^. 



(60) 



Tl = 



Moreover, we have 



1 / 1 \ , 



n=0 
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so that straightforward computation yields 



N 



N-l 



N 



(61) 



n=0 



So from (60 1 and (61) we deduce that 



N 



J I mK 

Aig^f(p= ■^,q= —jr> 



So clearly the discrete left-invariant vector fields acting on the discrete Gabor-transforms converge to 
the continuous vector fields acting on the continuous Gabor transforms pointwise as iV — > oo. 

In our algorithms it is essential that one works on the finite group with corresponding left-invariant 
vector fields. This is simply due to the fact that one computes finite Gabor-transforms (defined on 
the group f)r) to avoid sampling errors on the grid. Standard finite difference approximations of the 
continuous left-invariant vector fields do not appropriately deal with phase oscillations in the (discrete) 
Gabor transform. 

Remark 15 In the PDE-schemes which we will present in the next sections, such as for example the 
diffusion scheme in Section {Tl\ the solutions will leave the space of Gabor-transforms. In such cases one 
has to apply a left-invariant finite difference to a smooth function U G h2{Hr) defined on the Heisenberg- 
group Hr or, equivalently, one has to apply a finite difference to a smooth function U £ L2(M^) defined 
on phase space. If U is not the Gabor transform of an image it is usually inappropriate to use the final 
results in (57) and tSm on the group H,.. Instead one should just use 



{AY^U)[l,m,k] = (d7^^^[l,0,0]C/)[^,m,fc] 

{AY~U)[l,m,k] = (d7^^"[l,0,0]^7)[^,m,fc] 
{Afu)[l,m,k] 
iA^-U)[l,m,k] 



{dn^ [0,l,0]U)[l,m,k] 



(d7^^ [0, l,0]C/)[/,TO,fc] 



Ulll.7n.k]ll,0.0]]-U[l,m,k] 
[/[;,m,fc]-[/[[;,m,fc][- 1.0,0]] 
Ulll,m,k]lO,lfl]]-Ull,m,k] 
e/[;,m,fc]-^|^,m,fc] [0,-1,0]] 



(62) 



which does not require any interpolation between the discrete data iff ^ G N. The left-invariant operators 

'-). For example, {Af^U)[l,m] = [S^ o 



on phase space (58) and (59) are naturally extendable to 1^2{ 
o{S^y^U][l,m]^ K{e'^U[l + l,m] 



U[l, m]) for all U G £2({0, . . . , X - 1} x {0, . . . , M - 1}). 



9 Algorithm for the PDE-approach to Differential Reassign- 
ment 

Here we provide an explicit algorithm on the discrete Gabor transform G^i of the discrete signal f, 
that consistently corresponds to the theoretical PDE's on the continuous case as proposed in [B], i.e. 



convection equation (41 1 where we apply the first choice ([42|). Although that the PDE studied in [B] 



is not as simple as the second approach in (42) (which corresponds to a standard erosion step on the 
absolute value followed by a restoration of the phase afterwards) we do provide an explicit numerical 
scheme of this PDE, where we stay entirely in the discrete phase space. 

It should be stressed that taking straightforward central differences of the continuous differential 
operators of section [6] does not work, due to the fast oscillations (of the phase) of Gabor transforms. We 
need the lef-invariant differences on discrete Heisenberg groups discussed in the previous subsection. 

Explicit upwind scheme with left-invariant finite differences in pseudo-code for M = 1 

For I ^ 1, . . . , K - 1, m ^ 1, . . . M - 1 set W[l,m,0] := G:^f[/,m]. 
For i = 1,. . . ,T 

For « = 0, . . . , a: - 1, for m = 1, . . . , M - 1 set_^ 
v^[l,m]~~^ (log \W\[l l,m,t = 0] - log\W\[l - l,m,t = 0]) 
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«2 , m] — - ^ (log j IV ] , m + 1 , t = 0] 1 - log 1 1 [/ , m - 1 , f = 0] ) 

W[l,m,t] ■- W[l,m,t~l]+ K M [z+ {i^)[l,m] [A?' W][l,m,t\ + z- {v'-)[l,m] [A?^ W][l, 771,1]^ + 
MAt (z+{v^)[l,m] [A^'W][l,m,t\ + z-{i'^)[l,m] [A^^ W][l,m,t\ 



Explanation of all involved variables: 



I discrete position variable I = Q, . . . , K — 1. 

m discrete frequency variable m = 1, . . . , M — 1. 

t discrete time t = 1, . . .T , where T is the stopping time. 

ip discrete kernel ip = ij}^ = {V'a(riA^~^)}^ri(iv_i) or -0 = {i/'f N}^ri(^_i)See below. 

G^f[i,m] discrete Gabor transform computed by diagonalization via Zak transform [35] . 

W[l,m,t\ discrete evolving Gabor transform evaluated at position frequency m and time t. 

Af forward (+), backward (-) left-invariant position {i = 1) and frequency {i — 2) shifts. 

z^ z'^ ((f>)[l,m,t] — ma.x{(j)(l,m,t),0}, z~ {(j))[l,m,t\ = mm{(f>{l,m,t),0} for upwind. 



The discrete Cauchy Riemann kernel tp^ is derived in 7 and satisfies the system 

V,=o,...,K-iV™=o,....M-iVf6,,(,) : ^{A?^+Ar)+ia{A?^+A?')iG^or)[l,m]^0 , (63) 
which has a unique solution in case of extreme oversampling K = RI = N , L = 1. 

9.1 Evaluation of Reassignment 

We distinguished between two approaches to apply left-invariant adaptive convection on discrete Gabor- 



transforms (that we diagonalize by discrete Zak transform [35], recall subsection 7.1 1. Either we apply 



the numerical upwind PDE-scheme described in subsection [9] using the discrete left-invariant vector 



fields (58), or we apply erosion (44) on the modulus and restore the phase afterwards. Within each 



of the two approaches, we can use the discrete Cauchy-Riemann kernel or the sampled continuous 
Cauchy- Riemann kernel xp'^ . 

To evaluate these 4 methods we apply the reassignment scheme to the reassignment of a linear chirp 
that is multiplied by a modulated Gaussian and is sampled using N = 128 samples. A visualization of this 
complex valued signal can be found Fig. [4] (top) . The other signals in this figure are the reconstructions 
from the reassigned Gabor transforms that are given in Fig. [6] Here the topmost image shows the Gabor 
transform of the original signal. One can also find the reconstructions and reassigned Gabor transforms 
respectively using the four methods of reassignment. The parameters involved in generating these figures 
are N = 128, K = 128, M = 128, L = 1. Furthermore a = 1/6 and the time step for the PDE based 
method is set to At = 10~^. All images show a snapshot of the reassignment method stopped a,t t = 0.1. 
The signals are scaled such that their energy equals the energy of the input signal. This is needed to 
correct for the numerical diffusion the discretization scheme suffers from. Clearly the reassigned signals 
resemble the input signal quite well. The PDE scheme that uses the sampled continuous window shows 



Input signal 

e^.f(o = Rf (0+2^/(0 




Absolute value of 
Gabor transform 

\Qi,fll>,q)\ = \W^f{p,q,s)\ 



Absolute value of re-assigned 
Gabor transform 




p I. 



p <-)• I 



Figure 3: Illustration of reassignment by adaptive phase-invariant convection explained in Section [6j 
using the upwind scheme of subsection |9] applied on a Gabor transform. 
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ei 


£2 
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Erosion continuous window 


2.41 


10- 


-2 


8.38 


• 10- 


-3 


0.1 


Erosion discrete window 


8.25 


10- 


-2 


7.89 


• 10- 


-2 


0.1 


PDE continuous window 


2.16 


10- 


-2 


2.21 


• 10- 


-3 


0.1 


PDE discrete window 


1.47 


10- 


-2 


3.32 


• 10- 


-4 


0.1 




2.43 


10- 


-2 


6.43 


• 10- 


-3 


0.16 



Table 1: The first column shows ei = (|| f — f | 



12(1)) 



the relative error of the complex valued re- 



constructed signal compared to the input signal. In the second column £2 = (|| If] — |f| IU2(^))l|f|l7^ can 
be found which represents the relative error of the modulus of the signals. Parameters involved are 
K — M = N = 128, window scale a = | and convection time t = 0.1, with times step At = 10"'^ if 
applicable. PDE stand for the upwind scheme presented in subsection |9] and erosion means the morpho- 
logical erosion method given by eq. ( 44 ) . 



some defects. In contrast, the PDE scheme that uses i/j^ resembles the modulus of the original signal 



the most. Table pT] shows the relative £2-errors for all 4 experiments. Advantages of the erosion scheme 
( 44 1 over the PDE-scheme of Section |9] are: 

1. The erosion scheme does not produce numerical approximation-errors in the phase, which is evident 
since the phase is not used in the computations. 

2. The erosion scheme does not involve numerical diffusion as it does not suffer from finite step-sizes. 

3. The separable erosion scheme is much faster. 

The convection time in the erosion scheme is different than the convection time in the upwind-scheme, 
due to violation of the Cauchy-Riemann equations. Using a discrete window (that satisfies the discrete 
Cauchy-Riemann relations) in the PDE-scheme is more accurate. Thereby one can obtain more visual 
sharpening in the Gabor domain while obtaining similar relative errors in the signal domain. For example 
t = 0.16 for the PDE-scheme roughly corresponds to erosion schemes with t = 0.1 in the sense that the 
^2-errors nearly coincide, see Table 1. The method that uses a sampled version of the continuous window 
shows large errors and indeed in Fig.lHlthe defects are visible. This shows the importance of the window 
selection, i.e. in the PDE-schemes it is better to use window ■j/'^j rather than window ■j/'jj . However, 
Fig. [5] and Table 1 clearly indicate that in the erosion schemes it is better to choose window i/j^ than 

10 The Exact Analytic Solutions of Reassigned Gabor Trans- 
forms of ID-chirp Signals 

In the previous section we have introduced several numerical algorithms for differential reassignment. We 
compared them experimentally on the special case where the initial ID-signal (i.e. d = 1) is a chirp. In 
this section we will derive the analytic solution of this reassignment in the Gabor domain. Furthermore, 



we will show the geometrical meaning of the Cauchy-Riemann relation (381 in this example. 

Lemma 16 Let t > 0, rj G [|, 00). Let $^ : Tin C!{H{3)) be the operator that maps a Gabor transform 
W^^f to its reassigned Gabor transform 

Qt,a,^(|Wv,„/|)e"'-3{^*»/> 

where Qt.a.ri maps the absolute value |W^/,^/| of the Gabor transform to the unique viscosity solution 
Qt,ri{\yVipaf\) ■— W{-,t) of the Hamilton- J acobi equation 



'^{p.q,t) = -l^(a-^^{p,q,t)+a^^{p,q,t)) , (p, g) G i > 0, 
W{p,q,Q) = \yV^J{p,q,s)\ = \g^,{p,q)l {p,q) G 



(64) 



at fixed time t > 0. Set S^' :— o W^,„. Let Da be the unitary scaling operator given by [35). Then 



(-"'''/)(P,'7,'S) = Va(2^''X'„-i/)(a ^p,aq,s) (65) 
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10 20 30 40 50 60 70 80 90 100 110 120 



Figure 4: Reconstructions of the reassigned Gabor transforms of the original signal that is depicted on 
the top left whose absolute value of the Gabor transform is depicted on bottom left. In the right: 1st 
row corresponds to reassignment by the upwind scheme (A4 = 1) of subsection^ where again left we 

C D — 

used and right we used ip^^ . Parameters involved are grid constants K = M = N = 128, window 
scale a = 1/6, time step At = 10~^ and time t = 0.1. 2nd row to reassignment by morphological erosion 
where in the left we used kernel xp'^ and in the right we used xp^ . The goal of reassignment is achieved; 
all reconstructed signals are close to the original signal, whereas their corresponding Gabor transforms 
depicted in Fig. [6] are much sharper than the absolute value of the Gabor transform of the original 
depicted on the bottom left of this figure. 



Origiual 




Figure 5: The modulus of the signals in the bottom row of Fig. |4j For erosion (44) tp^ performs better 
than erosion applied on a Gabor transform constructed by ip^ ■ 




Figure 6: Absolute value of reassigned Gabor transforms of the signals depicted in the right of Fig. |4j 
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Proof First of all one has by substitution in integration that 

{W^J)ip,q,s) = Va(Wv,„.,I?a-i/)(^,a9,s) . (66) 

for ah / e L2(M'*), {p, q, s) G H{3). So if we introduce the scaling operator Da : L2(i/(3)) h2{H{s)) 
by 

(T>aU){p,q,s) ^ U{a ^p,qa,s),a>0, 



then (66) can be written as = y^DqW^j and by the chain-law for differentiation one has Qt,a,T] = 
DaQt, i.jjDa-i (where we note that the viscosity condition, cf. [26J is scaling covariant). Consequently, 
we obtain 

{E^-y){p,q,s) = ^e'-^-si^^afiP^^,^)} {BaQt,i,^Da-inaW^„V,-if){p,q,s) 
= ^e-g{w,,,/(a-ip,a,,s)} (Q,_,_^(2)^_i/))(a-ip,ag,s) 
= \/« (2^''/)(a" V, aq, s) □ 



Corollary 17 By means of the scaling relation Eq. {76) we may as well restrict ourselves to the case 
a = 1 for analytic solutions. 

Lemma 18 Let g(^) = e'^^'+^^+S with Re{w) < and Im{w) > 0. Then 

where the complex square root is taken using the usual branch-cut along the positive real axis. 

Proof By contour integration where the contour encloses a the two-sided section in the complex plane 
given by the intersection of the ball with radius R with < arg(z) < — iarg(— it;) (positively) and 
— TT < arg(2:) < — TT H — I arg(— w) (negatively) then one has by Cauchy's formula of integration and 
letting i? — >■ cx) that 
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/arg(z) = -i arg(-u)) \I-W\l\w\ 

from which the result follows. □ 
Theorem 19 Let r,b> 0. Let f he (a chirp signal) given by 

/(e)=e-|^e"'-«\ (67) 
Then its Gabor transform (with V'(0 — e^'^^ ) equals 



W^fip, 9, s) = J „ , , \ e-2-(^+¥)e(P.«)s.,.(p.9)" (68) 
V 27ro-^ — ?,r + 1 

with Brb = Re{Br^b) + iI'm{Br,b), Re{Br,b) < 0, {Im{BrM)V = I'miBr,b), iRe{Br,b)V = ReiBr,b) given 
by 

Orb ^254 + (_L+b2)2 I I ^^^4 _7r62(62_,_l 



-i(fe2 + J-) -Trr^b^ Trrfo^ 



Tirb^ + T + T^r^b^ 



and thereby we have for 77 € (5, 1]-' 

{'^t'''f)iP'1^s) ^ Vo'i^t^'^'^a-^f){a^^P,aq,s) ,a> 0, with 

(S^''/)(P,'Z,s) = e2-^'"^)e'MB..)(fc;'ee(-'-)««(^'-'')(-.-)^)(p,g) , ^ ' 
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with (positive) erosion kernel at time t > given by 

V - 
2r? 



knP,q) = + (70) 



for 77 S (5,1] a"-c? 

00 if p^ + > t^ 
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Proof Consider Eq. ([5|) where we set n = 1, d = 1, Eq. (67) and = er'"^ . We apply Lemma 
with w = iirr — —7t,u = —2'np — 2'Kiq, t — 2'Kipq — irp'^ and Eq. ( 68 ) foUows. Then we note that when 
transporting along equiphase planes, phase-covariance is the same as phase invariance and indeed by 
Lemma [l6| the main result (69) follows. Finally, wc note that the viscosity solutions of the erosion PDE 
is given by morphological convolution with the erosion kernel which by the Hopf-Lax formula equals 
k^{p,q) = t£^{t~^{p,q)'^) where the Lagrangian £jj{p,q) = {^H^){p,q) is obtained, cf. |26i ch:3.2.2], by 
the Fenchel-transform of the hamiltonian Hrj{a~^u,av) = ^(a^^(w)^ + 0^(1;)^)^'' that appears in the 
righthand side of the Hamilton- Jacobi equation Eq. (64), cf. [36, ch:2,p.24], Eq. (44), from which the 
result follows. □ 

Remark 20 In theorem we have set centered the chirp signal in position, frequency and phase. The 
general case follows by: {'WiM(po,qo.so)f){P^1^ = (Wv/)((Po, 9o, So)"^(p, <?, «))■ 

Remark 21 Note that trace{Re{Bj.b)} < and Det{Re{Brb)} > 0, so both eigenvalues of the symmetric 
matrix Re{Brb) are negative. Consequently, the equi-contours of the spectrogram (p,q) 1— >■ \G^f\{p,q) are 
ellipses. In a chirp signal ( without window, i.e. b ^ 00) frequency increases linear with ^ via rate r and 
thereby one expects the least amount of decay in the spectrogram along q = r p, i.e. one expects (1, r)^ to 
be the eigenvector with smallest absolute eigenvalue of Re{Br,b) ■ This is indeed only the case if b 00 
as 



Re{B, 



r.b) 



1 \ (52/2) // 1 



r 



r2fe4 + ((27r)-i+62)2 \\ r J ' \ 



1 



Remark 22 Note that Det{Re{Brb)} < 0, so the eigenvalues of the symmetric matrix Re(Brb) have 
different sign and consequently the equiphase lines in phase space (where we have set s — —pq/2) are 
hyperbolic. 

In case 77 — >■ 00 the Lagrangian is homogeneous, the Hamilton- Jacobi equation (that describes the 
evolution of geodesically equidistant surfaces {{p,q) G | W{p,q,t) = c}) becomes time-independent, 
cf. [Ml ch:4,p.l70] 

1 = a'^^W{p,q,t)+a^^W{p,q,t). 
dp'^ ap'^ 

In case = 5 the Hamiltonian is homogeneous and the erosion kernels are flat accordingly. The non- 
flatness of the erosi on k ernels can be controlled via parameter 77 € {^,00). Furthermore, the Hamilton- 
Jacobi equation in (64) is invariant under monotonic transformations iff 77 = ^ and this allows us to 

—1 - — 
compute {^^^ f ){p,q,s). 

Lemma 23 Let Ai,A2 denote the eigenvalues (with |Ai| < \X2\) of Re{Br.b) with respective normalized 
eigenvectors ki and Then each vector in can be written as (j>,q) — ai{p,q)ki + a2{p,q)k2. One 

has . 

Ai = ^ (-1 - 4b27r(l + 7r(&2 + r^)) + ^ (7rrb4)2 + 1 + f,2^(^2 _ ^,2)) = 

A2 = ^ (-1 - -ib^il + 7v{b' + r-2)) - ^{T:rb'iy + l(^J-+b2Tv{r2-b^)f ' 

with dbr = r^b"^ + (fo^ + and 

Cbr = v/1 + 8fe27r2(r2 + 62(-l + 27r2(64 + 262(2fe2 _ i)r2 + r^))) 
Nl = ^l+ '^-'^-^^rw^ ^^ for = 1, 2. 
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We omit the proof as the result follows by direct computation. 
Lemma 24 Consider the ellipsoidal isolines given by 

{{p',q') e I {p\q')BrAp',qT - Ai|aib',q')P + A2|«2(p',g')l' = c}, 
with c > arbitrary. Consider the line spanned by the principal direction with smallest absolute 



eigenvalue. Applying the settings of Lemma 23 this line is given by a2{p',q') = 0. Consider an- 



other point {p,q) on this line, i.e. a2ip,q) = p ^ rq, with p > 0. Then there exists a unique 
p* = p*{p,pt) = , ^ „ such that the circle {p — p*)"^ + (tp — rp*)^ — t^ is tangent to one of the ellipses 

{p' ,q')Br.b{p' iQ')^ — ^ilctiip' ,q')\^ + Xi\a2{p' , q')\'^ = c in {p*,Tp*). The corresponding time equals 



tji 



rXp.pr) = Vl + r2 ^/{p-p*{p,q)Y 
where the constants a and r are given by 



f^\p\ 



1 - 



a = b' 



et,^dt,^(-l + ct,r+46^7r^(6-r)(6+r))(l + e6^+4i)^7r(l + 7r(6^+r^))) 

((l + Cbr-4b47r2)2+8i)2-n-2(i + ci,^+4b<'(-l + 262),r2)r2 + 16b«'7r<'r'*)i (l-Ci,^+4ij2^(l+^(f,2_,_^2))) ' 



S|ll-l+4i,2^((,2_r2) • 



Proof In order to have the circle touching the ellipse in {p* ,Tp*), we have to solve the following system 

a2{p',q')=0, 

{p-p'f + {q-q'?=t\ 

^^{p'.q') = i,. 

The curvature along isocontours oiu{p,q) — {p' ,q')Br^b{p' ,q')^ is expressed as 



K{p',q' 



du(p' ,q') 8^u(p' ,q') _ 2 9u{p' ,q') dujp' ,q') d'^u(p' ,q') . / du{p' ,q') \ d'^ u{p' ,q' ) \ _ ^-L/Z/ , 
dq' c?p'2 Qq' Qp' dp' dq' \ dp' J dq''^ j 

then along ki we have a2{p',q') = i.e. q' = rp' and substitution in ( |72| ) yields 

^-'ip', rp') = (pT <y = t' = {p- pT + q' f = {P- + r') 



(72) 



this yields p'{p, rp) 



so since p* > p we find p* (p, pr 



□ 



of an eroded Gabor transform of a chirp 



The next theorem provides the exact solution Sj'^/ : 
signal at time t > 0. In the subsequent corollary we will show that the isocontours of the spectogram 
are non-ellipsoidal Jordan curves that shrink towards the principal eigenvector of Re(i?r,6) where they 
coUaps as t increases. This collapsing behavior can also be observed in the bottom rows of Figure [9] and 



Theorem 25 Let = | o,nd let t > 0. Let f be a chirp signal given by (61) then the reassigned Gabor 
transform of f (with window scale a= I) is given by 



27rb2-ir+l 
(At(p,9)) 



;(a2(P,<7) + t)^A2 



if ai{p,q)a2{p,q) ^ 0, 
ifai{p,q) = 0, 

if ct2{p, g) = and t < t^ax{p, rp). 



(73) 
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where Xi, X2 denote the negative eigenvalues (with |Ai| < \X2\) of Re{Br.b) and where {p, q) = ai{p, q)ki + 
a2{p,q)k2 with {fci, fc2} the normalized eigenvectors of Re{Br,b) and where the Euler- Lagrange multiplier 
Xt ■— Xtip, q) corresponds to the unique zero Pt{Xt) — of the A-th order polynomial 



4=1 



such that, for ai{p,q)a2{p,q) ^ 0, the corresponding unique solution {p' ,q') of 



{Re{Br,b) - Xt{p, q)) 



-^t{p,q) 



(74) 



(75) 



has maximum ^ Afc|afc(p', <?')p. The polynomial has at least 2 real-valued zeros. For each {p,q) £ 
fc=i 

there exists at — tmax{p, ?) ^ such that Pt has 3 real-valued zeros and Pg has 4 real-valued zeros for 
all s > t. Finally, the Lagrange multiplier satisfies the following scaling property 



Xt[p,q) = X^tiCpXq) for all C > 0. 



(76) 



Proof In case 77 = ^ we have 
{k!ef)ip,q) = 



mm 



^{p',q')RciB^_t)ip',q'y 



(77) 

Now Re(_Br.fc) < 0, hence the minimum of the associated continuous function can be found on the 
boundary of the convex domain. Apphcation of Euler-Lagrange yields the system 



\{p,q) 



p ~p 



(78) 



{p~p'Y + {q- q'Y = t 



First we must find the Euler-Lagrange multipliers Xt{p,q). If Xt '■= Xt{p,q) is not an eigenvalue of 
Ke{Br^b) the resolvent (Re{Br.b) — Xfl)^^ exists and one finds 



T\]2 



t' - ||Rc(S,,b)(Re(S,,b) - XtI)-\p,qY II 



which yields the 4th-order polynomial equation (74), where we note that {p,q) — Y^\^iCtk{p,q)^k with 
(ki , kj ) = 5ij so that 



t" = \\Rc{Br^b)iMBr,b) - Xtl)-\p,qf\f 



k=l 



>^kak{p, q) 



Xk - Xf {p,q) 



(79) 



Now that the (four, three or two) Lagrange multipliers are known, we select the one minimal associated 
to the global minimum and the general solution is given by 



,(At(p,9))2(p.g)(Rc(S,,,i,)-At/)-iRc(S,,b)(Ro(B,,6)-At/)-i 



(80) 



and substitution of {p,q) ~ J^j^i cgiP j q) ^i straightforwardly yields the first case in (73). The scaling 
relation (76) now directly follows by (80) and ai{(^p,(q) — (ai{p,q), for all i — 1,2, (p, 5) e and 

C>o. 

In case Xt{p, q) is equal to an eigenvalue of Re(i3r,6), say A; then (78 ) yields XkOtk{p' , q') = Xt{p, q)(ak(p' , q')- 



ak{p, q)), k — 1,2, from which we deduce that Eq. ( [74| ) is still valid if the resolvent does not exist, since 
then Xt{p, q) = Xi ^ ai{p, q) = 0. If ai{p, q) = then (p, q) is aligned with the (3 — l)-th principle axis of 
the ellipsoids {{p, g) e | {p, q)Re{Br^b)iP, qV = c}, c> and the minimum e(*+l"3-'l)^^3-' is obtained 
at {p',q') = {p,q)-\-tsgn{a3-i{p,q))k3^i = (a3_/(p, q) + isgn(a3_;)) k3_/, where only for / = 2 we get the 
extra condition t < ti^s.^(p,pT). See the second row of Figure |8j where in the third column t is chosen 
slightly larger than t-^sj^{p,pT). One can see that the minimum moves along the straight main principal 
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Figure 7: The Lagrange multiplier Xt{p,q) for {b,r,t) — (1, 1,0.1) used in exact S^'^/, cf. Figure 



10 



the corresponding approximations, in Eq. 
have set c = ^ ) . 



( 82 1 for respectively t > 



and t < t„ 



and 



(where we 



direction until {p* (p, rp), Tp*{p, rp)) is reached at time tniax(p,P'''), where the single minimum is cut into 
two minima that evolve to the side. 

For the / = 2-case a2{p, q) — (and t < t„ji^x{p,pT)) and the I = 1-case ai{p, q) — 0, see respectively 
first and second row of Figure [sj In these cases we find (by Eq. (79)) the Lagrange multiplier 



a; = ^ At = (1 + 



t 



3-1 



(81) 



which shows us that we indeed have a continuous transition of the solution across the principle directions 
in Eq. (73 1. Finally, we note that along the ray 3 fi (1 + ii){p* {p,pT),Tp* {p^pr)) there do not 
exist global minimizers of the optimization problem given by (77). □ 



Remark 26 We did not succeed in finding more tangible exact closed form expressions than the general 
tedious Gardano for mula for the Lagrange-multipliers Af(p, q) which are zeros of the fourth order poly- 
nomial given by Eq. \ H). For a plot of the graph of {p,q) i~> Xt{p,q), for {b,r) = 1), see Figurelli 
These plots do suggest a close approximation of the type 



{ ( 



Xtip,q) ~ < 



Al 1 + 



Aa 1 + 



Clip, 1)1 +\ci2(.P,l) 



2/ l°2(P>9)l^ + l°l(P,<;)|-' \ 2 



for t,nax(p, q) <t 

for t 

max {p, q)>t , 



(82) 



with c G (0, 1], where we recall that tmax{p,Q) is defined in Theorem 25 This approximation obeys the 



scaling property (76) and is exact along the principal axes (i.e. exact if ai =0 or a2 — 0). 



Corollary 27 Let Ai,A2 < be the eigenvalues of Re{Br^b) with |Ai| < IA2I and with corresponding 

eigenvectors ki, k^. The function ^t'^f '. H,. — > C converges pointwise towards as i — >■ 00. The 

1 

isocontours of f \ with t > are non-ellipsoidal Jordan curves that retain the reflectional symmetry 
in the principal axes of Re{Br,b)- The anisotropy of these Jordan curves (i.e. the aspect ratio of the 
intersection with the principal axes) equals 

(83) 

- 1 

which tends to 00 ast"[tfin ■= which is the finite final time where the isocontour {{p,q) \ |^('^/(p, g, s)| 

c} collapses to the span of k 



Proof Set rj = ^. By the semi-group property of the erosion with kernel Eq. (71 1 and Eq. (44) (i.e. 
~ ^t+s) ^'^'i the fact that the function |Hq' ^ f\ = \G^f\ vanishes at infinity and has a single critical 
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Figure 8: Various plots of the solutions of optimization problem j??! arising in Sj'^. The horizontal 
axis, is the p' axis, whereas the vertical axis is the g'-axis. Top row (p, q) — (0.02, 0.02) — 0.02ki, middle 
row (p, q) = (0.05, —0.05) = — 0.05k2, bottom row {p, q) — (—0.1, —0.05) for increasing time-frames t > 
from left to right. The minimizer(s) on the circle is(are) indicated in Green, whereas the other spurious 
solution(s) of Eq. (75) and Eq. (74) is(arc) indicated in black. 



point; a maximum at (p, q) — (0, 0) it follows that 

I'^'t+s/l = kg Q I'E.l'^ f\ < attains a single extremum; that is a maximum at {p,q) = (0,0), 

|Sj'''/| is continuous , |Sj'''/|(p, (?) ^ as \\{p,q)\\ oo, 

for all t,S > 0. Consequently, the isocontour of with value c > is a Jordan curve that is strictly 

contained in the interior of the isocontour of with the same value. By Theorem 
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the isocontours 

are ellipsoidal for t ~ 0. For t > these are no longer ellipsoidal since {p, q) i— >■ At(p, q) is not constant in 



Eq. (73). For i — > oo the kernel fc^ , cf. Eq. (71), vanishes and the erosion, cf. Eq. (44), with the kernel 



k^ converges for each point (p, q) G to the global minimum of \G4>f \ which is zero. 

Eq. ( 74 1 is invariant under ai i— > —a , and thereby At is invariant under reflections in the principal 
axes of Re(i3r,&)- As a result, by Eq. (73), all isocontours of are invariant under these reflections. 

Such an isocontour is given by 



\ai{p, g)pAi 



1^2 (p, q)?\2 



(Ai-At(p,g))2 (A2-At(p,<7))^ 



(84) 



for some C < 0. Now set ai = in Eq.(84) and solve for a2, then set a2 = in Eq.(84| and solve for 
ai. Here one should use the exact formula, Eq. (81) for A* that holds along the principal axes. Finally, 
division of the 2 obtained results yields Eq. (83 1 from which the result follows. □ 



See Figure [s] for solutions of the optimization problem (77 1 for various settings of {p,q,t) 



10.1 Plots of the Exact Re-assigned Gabor Transforms of Chirp Signals 



For a plot of the exact ^t^f in Theorem 
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usmg Gardano's formula for the zeros of the fourth order 
polynomial given by Eq. (74) we refer to Figure^ (for the case (&,r) = (^,1)) and Figure 



10 



(for the 

case (6, r) — (1, 1)). Withinthcse Figures one can clearly see that during the erosion process isocontours 
collapse towards the eigenspace (ki) with smallest eigenvalue. Moreover for (fe, r) = (1,1) one has 
ki « ^)"^ s-nd finally, we note that the discrete Gabor transforms closely approximate the exact 



Gabor transforms as can be seen in the top row of Figure 10 Although that, for 77 = ^, the eroded 
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chirp 



discrete GT 




K = N = M = 128 
a = L=l 



(r,6) = (l,l) 5 ^ 128 

phase continuous GT modulus continuous GT evolving over time 




Figure 9: Visualization of {p, q, s) H' (^J'^ g, s) where / is the chirp signal (67) with (6, r) = 1). 



Top row: left; the input chirp signal, right: We have sampled the chirp signal onT28-grid and depicted 
the corresponding discrete Gabor transform {K = M = N,L — a — 1) where color represents phase and 
where intensity represents the modulus. Middle row; equisurfaces of the phase and modulus of our exact 



solutions for =1'^ f in Theorem 25 phase is constant over time whereas the modulus is considered for 
t = 0, < = 0.1 and t = 0.2. Iso-intensities have been fixed (-0.05, -0.09, -0.14, -0.18, -0.23) over time 
to show their collapsing behavior towards ki in accordance with Corollary [27j Bottom row; Restrictions 
to the surface s = — ^ yields the solutions in phase space. 



ellipses are no longer ellipses the principal axis (ki) is preserved as a symmetry- axis for the iso-lines in 
phase space (where s — — ^). 

11 Left Invariant Diffusion on Gabor transforms 

A common technique in imaging to enhance images / : — >■ K via non-linear adaptive diffusions are 
so-called coherence enhancing diffusion (CED) schemes, cf. [371 [53 El] > where one considers a diffusion 
equation, where the diffusivity/conductivity matrix in between the divergence and gradient is adapted 
to the Gaussian gradient (structure tensor) or Hessian tensor computed from the data. Here the aim 
is to diffuse strongly along edges/lines and little orthogonal to them. In case of edge-adaptation the 
anisotropic diffusivity matrix is diagonalized along the eigenvectors of auxiliary matrix 

A{uf{;s)){x,y) ^ {G,*iVufi;s){Vuf{;s)f))ix,y) (85) 

where Vu / denotes either the Gaussian gradient of image / or the gradient of the non-linearly evolved 

image m/(-,s) and where the final convolution, with Gaussian kernel G^{x,y) = (47re)~^e~T^ with 
e > small, is applied componentwise. In case of line-adaptation the auxiliary matrix is obtained by the 
Hessian: 

(A(u/(-,s)))(x) = (G,*i7u/(-,s))(x) = [(Ge*a^.a^^.w/(-,s))(x)],,j-=i^2 = [{d^^d^^G,*uf{-, s)){ji)]^^j=i,2 ■ 

(86) 
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Given the eigen system of the auxiliary matrix, standard coherence enhancing diffusion equations on 
images cf. [37] can be formulated as 





(1 - e)e "^-^2 



^ ^ V (1 - £)e (^1-^2)^ +e y V y 

u/(x, 0) = /(x), xeB^. 

Here we expressed the diffusion equations in both the global standard basis {e^;, e^} ;= {(1, 0), (0, 1)} 



{dx,dy} and in the locally adapted basis of eigenvectors of auxiliary matrix A{uf{-, s))(x), recall Eq. (85 ) 



and Eq. (86) 



{ei,e2} {ei (A(w/(-, s))(x)) , ea (A(w/(-, s))(x))} 4-^ {8^,0^}, 

with respective eigenvalues A^- :— Xk{A{uf{-,s)){x)), k = 1,2. The corresponding orthogonal basis 
transform which maps the standard basis vectors to the eigenvectors {61,62} is denoted by S = (ei | 62) 
and we have {du dy) — {dx dy) ■ S. Note that 

A = S-diag{Ai,A2}-S-i. 

At isotropic areas Ai — >■ A2 and thereby the conductivity matrix becomes a multiple of the identity 
yielding isotropic diffusion only at isotropic areas, which is desirable for noise- removal. 

A typical drawback of these coherence enhancing diffusion directly applied to images, is that the 
direction of the image gradient (or Hessian eigenvectors) is ill-posed in the vicinity of crossings and 
junctions. Therefore, in their recent works on orientation scores, cf. [40l [121 [10] Franken and Duits 
developed coherence enhancing diffusion via invertible orientation scores (CEDOS) where crossings are 
generically disentangled allowing crossing preserving diffusion, see Figure |11| The key idea here is to 
extend the image domain to a larger Lie group SE{2) = x SO (2) where the left-invariant vector 
fields provide per spatial position a whole family of oriented reference frames cf. [HI |H] . This allows the 
inclusion of coherent alignment ( "context" ) of disentangled local line fragments visible in the orientation 
score. 

In this section we would like to extend this general idea to Gabor transforms defined on the the 
Heisenberg group, where the left-invariant vector fields provide per spatial position a whole family of 
reference frames w.r.t. frequency and phase. Again we would like to achieve coherent alignment of all 
Gabor atoms in the Gabor domain via left-invariant diffusion. 
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original image 



CED 



CEDOS 




Figure 11: Left-invariant processing via invertible orientation scores is the right approach to generically 
deal with crossings and bifurcations. Left column: original images. Middle column: result of coherence 
enhancing diffusion on images (CED), cf. 1411 . Right column: coherence enhancing diffusion via orienta- 
tion score, cf. [12l|40]. Top row: 2-photon microscopy image of bone tissue. 2nd row : collagen fibers of 
the heart. 3rd row: artificial pattern, coherence enhancing diffusion on orientation scores (CEDOS) is 
capable of handling crossings and bifurcations, whereas (CED) produces spurious artifacts. 



In order to generalize the CED (coherence enhancing diffusion) schemes to Gabor transforms we must 
replace the left-invariant vector fields {dx, dy} on the additive group (IR^, -|-), by the left-invariant vector 
fields on Hr- Furthermore, we replace the 2D-image by the Gabor transform of a ID-signal, the group 
by Hr, the standard inner product on by the first fundamental form (37) parameterized by (3. Finally, 



we replace the basis of normalized left-invariant vector fields {e^., ej,} := {(1, 0), (0, 1)} O {dx,dy} on 
by the normalized left-invariant vector fields {e^;, e^} :— {(1, 0), (0, 1)} O ^2} on Hr. 

These steps produce the following system for adaptive left-invariant diffusion on Gabor transforms: 



dtW{p,q,s,t) =i(i3-^-Ai A2)-S- 



(l-e)e- 



13-^ AlW 

A2W 



(P, 9, s, t) 



~ " " > (l-e)e (^1-^2)^ + 
for all (p, q, s) £ H^, t > and fixed c > 0, e > 



W{p, q, s, t). 



W(p, q, s, 0) = W^/(p, q, s), for all (p, q, s) £ Hr. 

with S := (ei I ea) and {duW dyW) = il3-^AiW A2W)S and where 

Afe :=Afc(A(|W^/|)(p,Q,s)),fc = l,2,|Ai| < IA2I 
efc := efc(A(|W^/|)(p,g,s)),fc = l,2, 



(89) 



denote the eigenvalues and {61,82} O {du,dy} the corresponding normalized eigenvectors of a local 
auxiliary matrix A(|yV^/|)(p, q, s): 



A(|>V^/|)(p,g,s) ek = Xk efc. 
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The eigenvectors are normalized w.r.t. first fundamental form (37) parameterized by /3, i.e 



efei = e^e2 = e^e^ — e^ey = 1 => = S^^ 

This auxiliary matrix A(|W^/|)(p, q, s) at each position {p,q,s) G Hj. is chosen such that it depends 
only on the absolute value |W^/| (so it is independent of phase parameter s) of the initial condition. 
The same holds for the corresponding conductivity matrix- valued function appearing in Eq. (88): 



e 
(1 - e)e~'^I^ +£ 

where we note that A = S • diag{Ai, A2} • S. 

We propose the following specific choices of auxiliary matrices 

A{\W^f\){p,q,3) = [{dp^dp^G, * mf\){p,q)kj=l.2 , 
Ai\W^f\){p,q,s) = [{dp^G^ * \W^mp,q) idp^G,*\W^f\){p,q)],.^j^i^2 



with pi 



^p and P2 — q and separable Gaussian kernel G^ {p, q) 



Note that 



Ai\W^,f\_^ dp\W^,f \ and A2\Wy,f \ = dg\W^f\. Now since |W^/| = \g^f \ and A = S o A, o S'^ 
(recall (29)) and S o W^f = G^pf we get the following equivalent non-linear left-invariant diffusion 
equations on phase space: 



dtW{p,q..t) 



W{p,q,0) 



( r^Ai A2 ) • s- 
( a,, a, ) 



e 

(l-e)e"^ 

£ 

(1 - e)e '^1-^2)^ + e 
for all (p, q) e R^, t > and fixed c > 0, £ > 0, 

- G^fip,q): for all e R". 



A2W 



{p,g,t) 



W{p.q..t), 



(90) 



with again S = (ei | 62) and {duW dyW) = (/3 A\W A2W) S and where we recaU from (89) that Afc, 



efc denote the eigenvalues and normalized eigenvectors of the local auxiliary-matrix A(|yV^/|)(x, s) 
A(|5^/|)(x, s). In phase space, these eigenvectors (in M^) correspond to 

{ei,e2} o {du,dij} = {5 o a„ o S^^,S o dy o S^^} , 



so that we indeed get the right correspondence between (90) and (88): 



V(p,g,s)effr,t>0 
V(p,g,s)eff,,,t>0 
V(p,g,3)eff^,t>0 

^{p,9)eR2,t>0 ■ 



dtW{p,q,s,t) = ( 



(1 - £)e (^1-^2)^ + 
f e 

dt(s-^w)(p.q,s.t) = { a„ a„ ) ^ ^ 

^ ^ ^ ^ ^ 1^ (1 - £)e '^1-^2)^ + 



) 



(1 - £)e (^1-^2)^ + 



(1 - £)e (Ai-A2)^ ^ J 

J W{p,q,t) 



j W(p,9,s,t) 

5 ( (5-il4')(p,<7,s,t) « 



a,. 



so that the uniqueness of the solutions W oi ( 90 ) and of ( 88 ) implies that 



^^(•,•,0) = g^/ = 5oWv,/ = 5oW^(.,.,.,0)=^Vt>o : W{-,-,t) ^ S oW{-,-,-,t) 



Clearly, Problem (90) is preferable over problem (88) if it comes to numerical schemes as it is a 2D- 
e volution. 



See Figure 12 for an explicit example, where we applied adaptive nonlinear diffusion on the Gabor 
transform of a noisy chirp signal. Due to left-invariance of our diffusions the phase is treated appropriately 
and thereby the effective corresponding denoising operator on the signal is sensible. Moreover, the tails 
of the enhanced chirp are smoothly dampened. 
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Figure 12: Illustration of the left-invariant diffusions on Gabor transforms that are adapted to the 
Hessian of the absolute value. The corresponding operator / i-> Q*^^j^Q^, in the signal domain smoothly 
enhances the signal without trcsholding of Gabor coefficients. In the most right plot we depicted ellipsoids 
representing the local eigenvectors of the Hessian of along which the diffusion locally takes place. 
The directions of the ellipsoids coincide with the eigenvectors of the Hessian, whereas the anisotropy of 
the ellipsoids is determined by the fraction IA1/A2I of the eigenvalues {Ai, A2} with |Ai| > IA2I, the colors 
indicate the directions of the largest eigenvector and the intensity reflects the relative damping factor 
{p,q) !->• 1 + e-i(l - e)e-'=(^i(P'«)-^2(p,9))"";s^^(p^q) of the second eigenvector. 
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Figure 13: Application of re-assignment of Gabor transforms of 2D-images. The dot indicates the 
position (p^jP^) e M.'^ in the image (left) where we depicted the local frequency distribution {q^^q^) i— )■ 
9^) (right) where we color-coded the phase. Top row input at two locations. Middle 
row output for t = 0.1. Bottom row output for t = 4. Further parameter settings: a ~ 1, N — 64, 
K = M ^ 32, L = 2, = D"^^ = 1. 

12 Reassignment, Texture Enhancement and Frequency Esti- 
mation in 2D-images {d = 2) 

In Section [6] we applied differential reassignment to ID-signals. This technique can also be applied 
to Gabor transform of 2D-images. Here we choose for the second option in ( [42| ) as the correspond- 
ing algorithm is faster. In this case the reassigned Gabor transforms concentrate towards the lines 
{{p\p\quq2,s) G i?5 I dp,\g^f\{p,q,s) = 5,. (p, g, s) = 0,z = l,2,p = {p\p'),q = (91,92) € K^}, 
which coincides with the stationary solutions of the corresponding Hamilton Jacobi equation on the mod- 
ulus. See Figure [TSj where the amount of large Gabor coefficients is strongly reduced while maintaining 
the original image. 

12.1 Left-invariant Diffusion on Phase Space as a Pre-processing Step for 
Differential Reassignment 

Gabor transforms of noisy medical images often require smoothing as pre-processing before differential 
reassignment can be applied. Linear left-invariant diffusion is a good choice for such a smoothing as a 
pre-processing step for differential reassignment and/or local frequency extraction, since such smoothing 
does not affect the original signal (up to a scalar multiplication). Here we aim for pre-processing before 
reassignment rather than signal enhancement. The Dunford-Pettis Theorem |42j shows minor conditions 
for a linear operator on L2(X), with X a measurable space, to be a kernel operator. Furthermore, all 
linear left-invariant kernel operators on L2(i?2d+i) are convolution operators. Next we classify all left- 
invariant kernel operators on phase space (i.e. all operators that commute with £(p ^ s) :— SC(^p_q^s)S^^). 

Lemma 28 A kernel operator ICk : L2(M2'^) hoo{R^'^) n L2(M2'^) given by 

{ICkU){p,q)= I l{p,q,p\c^)U{p',q'WAq\ 

with /c e Li ( left-invariant if 

k{p, q,p',q')= e^^^^(f ) fc(p + p,q + q,p' +p,q' + q) (91) 
for almost every (p, q), {p' , q'), {p, q) E M^'*. 
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Proof We have LgU{p' , q') = {SCgS-W){p', q') = e2^^''+'^*(29P'-P'3)f7(p' ~p,q' ~ q) so that 
(JCkCgtlM q) = /g,, ~k{p, q,p', g')e2--+"(2gp'-P9);7(p' -p,q'- q)dp'dq' 

from which the result follows by substitution p" = p' — p, q" = q' — q. □ 
Theorem 29 Let d G N. Let $t 6e a left-invariant semi-group operator on Hr given by 

($,(C/))(.9) = (Kt * U){g) / i^t(/i-\g)C/(/i) d^C/i) 



with fx the left-invariant measure on Hr, t > 0. Then the corresponding operator on phase space = 
5 o $4 o is given by 

{^t{g^f)){p,q) = / / kt{p,q,p\q')g^fip',q')dp'dq\ 

(92) 

with kt{p, q,p', q') = / e-"^P'-9'X,(p _ g _ q', - s' _ i(p . q'-q . p'))e~^-'^^' ds'. 



d 2d 

In particular we consider the diffusion kernel for horizontal diffusion generated by D^^ ^ (yli)2+I?22 ^ (^i)^ 

i=l i=d+l 

on the sub-Riemannian manifold {H2d+i, ds + • dg — g • dp)). 



Proof Eq. (92) follows by direct computation where we note h ~ {p',q',s') = {—p',—q',~s') and 
taking the phase inside, which is possible since {SW^f){p' ,q') — {G4,f){p',q') is independent of s' . The 
heat kernel on Hr is obtained by the heat kernel on H2d+i via K^'' {p,q, s) — ^ K^^{p,q,s -\- k). □ 

As the analytic closed form solution of this kernel is intangible, cf. [2H1IIH]) we consider the well-known 
local analytic approximation 

r 1 .('^^^^ , IIpIIVII'II^^ 

Kt{p,q,s) = , -f^ l,v'IJTTij22+-5Tr + 3^;4t_ ^gg^ 

StVD^^D^^ (47rte-iDii)i(47rte-iL»22)f 

This approximation is due to the ball-box theorem |43j or theory of weighted sub-coercive operators on 
Lie groups [44J, where we assign weights Wi = 1, i = 1 . . . ,2d and W2d-\-i = 2 to the left-invariant vector 
fields {^1, . . . ,^2<i+i}- 



Lemma 30 Let Kt : H{2d + 1) — > be given by Eq. (93). Then the corresponding kernel on phase 
space is given by 



Kt[P,q,P,q) e l + 647?iiD22t2c-2- ^y*^ 



and satisfies Eq. (91). 



Proof By Eq. (93) and Eq. (92) we find by substitution v = s' — s — k + \{p ■ q' — q ■ p') that 

l,t[P,q,P,q) e e 8t^D^ ^i,t,-,D^,^^^i,tc-^o--}i 

l+s'-fc+i(p-g'-(j-p') _|^|_ 



The integral can be computed by means Jj^e '""''''e dt; = ^2 ^4^2 for 7 > . The kernel in Eq.(94) 
satisfies the left-invariance constraint ([9T|) since —{p~p){q + q) = 2q{p — p) — (p — p){q + q + 2q) 13 
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Figure 14: Top row: Input image and two slices |^^/i(Pfc, •, ■){, k — 1, 2, in the Gabor domain centered 
around the indicated positions. Bottom row: output in spatial and Gabor domain of linear diffusion 
combined with a squaring of the modulus. (Left: input image /i, i.e. i = 1, right: input image /2, 
i.e. i — 2). Akin to left-invariant diffusions on orientation scores where we can generically diffuse along 
crossing lines (Recall Figure 111 we can apply left-invariant diffusion on Gabor transforms and thereby 
diffuse along crossing frequencies. 



12.1.1 Possible Extension to Texture Enhancement in 2D images via Left-invariant Evo- 
lutions 

The techniques signal enhancement by non-linear left-invariant diffusion on Gabor transforms in subsec- 
tion 11 can be extended to enhancement of local 2D-frequency patterns and/or textures. This can have 
similar applications as the enhancement of lines via non-linear diffusion on invertible orientation scores, 
cf. (12] and Figure [IT] However, such an extension would yield a technical and slow algorithm. In- 
stead, akin to our earlier works on contour enhancement via invertible orientation scores, cf. [14| , lll j . we 
can use the (more basic) concatenation of linear left-invariant diffusion and monotonic transformations 
in the co-domain. Figure [Ml shows a basic experiment of such an approach. 



12.2 Local Frequency Estimation in Cardiac Tagged MRI Images 

In the limiting case differential reassignment concentrates around local maxima of the absolute value of a 
Gabor transform. These local maxima produce per position {p^,p^) e a local frequency {qi, (72) G 
estimation. The reliability of such local frequency estimations can be checked via differential reas- 
signment. In this section we briefly explain a 2D medical imaging application where local frequency 
estimation is important for measuring heart wall deformations during a systole. For more details, see 

m- 

Quantification of cardiac wall motion may help in (early) diagnosis of cardiac abnormalities such 
as ischemia and myocardial infarction. To characterize the dynamic behavior of the cardiac muscle, 
non-invasive acquisition techniques such as MRI tagging can be applied. This allows to locally imprint 
brightness patterns in the muscle, which deform accordingly and allow detailed assessment of myocardial 
motion. Several optical flow techniques have been considered in this application. However, as the 
constant brightness assumption [JSJ 271 155^ in this application is invalid, these techniques end up in a 
concatenation of technical procedures, |70, 49, 50 yielding a complicated overall model and algorithm. 
For example, instead of tracking constant brightness one can track local extrema in scale space |51) 
(taking into account covariant derivatives and Helmholtz decomposition of the flow field f20^) or one can 
compute optical flow fields [501 123] that follow equiphase lines, where the phase is computed by Gabor 
filtering techniques known as the harmonic phase method (HARP), cf. [52] • Gabor filtering techniques 
were also used in a recent applied approach "1^] where one obtains cardiac wall deformations directly 
from local scalar-valued frequency estimations in tagging directions. Here we aim for a similar short-cut, 
but in contrast to the approach in |19) we extract the maxima from all Gabor transform coefficients 
producing per tagging direction an accurate frequency covector field / — fidx^ + f2'ix^ (not necessarily 
aligned with the tagging direction). From these covector fields one can deduce the required deformation 
gradient D by duality as we briefly explain next. 
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Let /j ' : — > M be the images corresponding to tag-direction 6i G [0, 27r), i — 1, . . . , N , with N > 2 
the total number of tags at time t > corresponding to (independent) tag-direction (cos(6'i), sin(6'i)) 
with 0i e [0, 27r). We compute the Gabor transforms G^{ft^) : M"' — >■ C and apply a linear left-invariant 
evolution and extract per position p the remaining maxima w.r.t. the frequency variable q. For each 
position p = (p^yP^) € the Gabor transform 5^(/f')(p, •) typically shows only two dominant and 



noisy blobs (for t small), cf. Figure 15 The 2 blobs relate to each other by reflection q i— )■ — q. 



Best frequency estimates are obtained by extracting maxima after applying a large time evolution. 
By the results in |53j this maximum converges to the center of mass, that gives a sub-pixel accurate 
estimate. This yields per position (p^,p^) € and per tagging direction 6i the local frequency estimate 
q'''(p^,p^) :— ±(g^'*(p^,p^), 92 Xp^jP^)) that we store in a covector field 

where da;^ = dx, dx^ — dy. The obtained frequency field q*'* can be related to the requested deformation 
gradient Dt = [ ^ ] , where xt denotes the position of a material point in the heart wall at time t > 0. 
We assume duality between frequencies and velocities by imposing 



(q*'^(xO, -^Ms) 
as 



(Xt-l), — Xt_i(s) 

as 



(95) 



for all smooth parameterizations (s,t) i— > xt(s) e with xt(0) ~ x^ and alH = 1, . . . , A^. This yields 



E g*^Xx0^i(xt-i) iti(O) = E (xO i^(0) = E q^-' 
Qt(xt)Dt(xt_i)xt_i = Qj_i(xt_i)xt^i 



'(x,_oiri(o)^ 



with Df = [Dl] = 



dxl 



9, 



e M^^2 a^^(j xt(0) (if ,a;?) = ^Xt(O) arbitrary so that 



Q,(xi)Dt(xt„i) = Q,„i(xf_i). 

Provided that the frequency estimates are smooth and slowly varying with respect to the position variable 
one can apply first order Taylor expansion on the left-hand side and evaluate Qj(x4_i) in stead of Q((xt). 
This corresponds to linear deformation theory, where Dt(xt) k, Dt(xt_i) provided that displacements 
||xt — xt_i|| are small. Consequently, we use per (fixed) spatial position x G K^, the least squares solution 
of Qt(x)Dt(x) = Qt_i(x) given by 



D,(x) = ((Q,(x))^Q,(x))-i(Q,(x))^Q,_i(x) 



(96) 



as our deformation gradient estimate. For experiments on a phantom with both considerable non-uniform 



scaling, fading, and non-uniform rotation see Figure 16 The iteratively computed deformation net is 
close to the ground truth deformation net. 

The deformation nets are computed as follows. We fix a single material point at the boundary (the 
green point in Figure 16 ) and first compute the grid points on the outer contour in circumferential 



direction. From these points we compute the inner grid-points iteratively in inward radial direction. Let 
Xj'-' denote the j-th grid-point in the deformation net at time t > on the r-th closed contour. Let 
J denote the number of points on one closed circumferential contour and let R denote the number of 
points on a circumferential contour, let T be the total time (number of tags), then 

Single material point x^'"^ given for alH e {0, ... T — 1}, 
Initial polar grid Xq-' given for all r e {1, . . . R},j G {1, . . . J}, 

Deformation field Dt computed for alH > from the local frequency fields {q*'*}^i via Eq. (961, 
Then for t = 1 to T-1 do 



for J = 1 to J compute x 
for r = 1 to R compute x^"''^'"' 



X, 



(x^^^'-x^f,) 



(97) 

In contrast to well-established optical flow methods, [351 [Ml [13 [SDl UHl 123] this method is robust with 
respect to inevitable fading of the tag-lines, and can be applied directly on the original MRI-tagging 
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Figure 15: Top row: MRI-tagging images (on a 64 x 64-grid) at time frame t = 2 and time frame t = 8 
(the systolic phase of the heart is sampled in 30 time frames), scanned in 4 tag-directions, i.e. = 4. 
Left: horizontal tag 6*1 = at i = 2 with estimated frequency field q^'^, right: vertical tag 6*3 = 7r/2 
with estimated frequency field q^''^ plotted on top. In the frequency estimates we used maximum over- 
sampling L = 1, = 64 in the discrete Gabor transform Eq. (471 with Gaussian window-size with a 
standard deviation of u = 6 (pixel lengths). Middle row: Corresponding deformation nets computed from 
the frequency fields and outer boundary via Eq. (96) and Eq. (97) Bottom row: Plots of the absolute 
value of the Gabor transform (91,92) \Gipf{p^^,qi,q2)\ of MRI-tagging image / with orientation 
02 — tt/4: restricted to the green and blue point {p^,p^) of interest. We have masked the low frequencies 
to show the 2 frequency blobs of the tagging lines. 
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Figure 16: Ground truth phantom where the middle image at time t > is obtained from the left image 
t = via both non-uniform scaling and non-uniform rotation. The phantom includes a strong fading 
of the tag-line intensities in this phantom. The deformation net depicted on top in red is computed 
according to Eq. (97), with N ^ 4, R = 8 and J = 50. The deformation net is entirely computed from 
a single indicated material point at the boundary. The right image shows a comparison between the 
computed deformation net (red) with the ground-truth deformation net (grey). 



images without harmonic phase or sine phase pre-processing [52l[50]. Some of the optic flow methods, 
such as [S5J [Sni [S7] , do allow direct computation on the original MRI-tagging images as well, but they 
are relatively expensive, complicated and technical, e.g. |20j . 

First experiments show that our method can handle relatively large nonlinear deformations in both 
radial and circumferential direction. The frequency and deformation fields look very promising from the 
qualitative viewpoint, e.g. Figure [16] and Figure |15[ Quantitative comparisons to other methods such 
as [ini m] and investigation of the reliability of the frequency estimates (and their connection to 
left-invariant evolutions on Gabor transforms) are interesting topics for future work. 

Although Gabor transforms are widely used in ID-signal processing, they are less common in 2D- 
imaging. However, the results in this section indicate applications of both Gabor transforms and the 
left-invariant evolutions acting on them in 2D-image processing. 

Acknowledgements The authors gratefully acknowledge Jos Westenberg (Leiden University Medical 
Center) for providing us the MRI-tagging data sets that were acquired in 4 different tagging directions. 



A Existence and Uniqueness of the Evolution Solutions 



The convection diffusion systems ( |17[ ) have unique solutions, since the coefficients and Dij depend 
smoothly on the modulus of the initial condition |yV^/| — \G^pf\- So for a given initial condition W^/ 
the left-invariant convection diffusion generator Q{\yVt(,f\,Ai, . . . ,A2d) is of the type 

d d 

Q{\W,pflAi, . . .,A2d) = ^a^A + ^AjAj. 

i—l — l 

Such hypo-elliptic operators with almost everywhere smooth coefficients given by ai(p, q) = ai{\Qipf\){p, q) 
and I3ij{p,q) = Dij{\Q^f\){p,q) generate strongly continuous, semigroups on L2(M^), as long as we keep 
the functions ai and Pij fixed, [581 130] . yielding unique solutions where at least formally we may write 



W{p,q,S,t) = <^>t{W^,f){p,q,s) = e\^-' JW^f{p,q,s) 

with limt^f) W {■ , t) — W^f in L2-sense. Note that if -0 is a Gaussian kernel and / 7^ the Gabor 
transform Q^f ^ is real analytic on M^'', so it can not vanish on a set with positive measure, so 
that ai : — !■ M are almost everywhere smooth. This applies in particular to the first reassignment 
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approach in (42 ) (mapping everything consistently into phase space using Theorem |4]), where we have set 
C = 0, aii\g^,n) = Mi\g^fm4J\-'dj,\g^f\ and a2{\g^f\) = Mi\g^,fm^f\-^d,\g^f\. in the second 
approach in (42 1 the operator U i— ?> C{U) is non-Unear, left-invariant and maps the space L^(IR^) = 
{/ e L2(K^) I / > 0} into itself again. In these cases the erosion solutions ( [44| are the unique viscosity 
solutions, of (41), see [59] . 



Remark 31 For the diffusion case, ^ ch:7], ^ ch:6j, we have D — [-Dy]ij=i,...,2!i > 0, in which case 
the (horizontal) diffusion generator Qd W^/|, ^i, . . . , y^2d) on the group is hypo- elliptic, whereas the 
corresponding generator Q(\g^f\,Ai, ... ,A2d) on phase space is elliptic. By the results f26, ch:7.1.1] 
and IMl there exists a unique weak solution W = SW G L2(M+,HIi(M2))nHIi(M+,L2(M2)). By means of 
the combination of Theorem^and Theorem^we can transfer the existence and uniqueness result for the 
elliptic diffusions on phase space to the existence and uniqueness result for the hypo-elliptic diffusions on 
the group (that is via conjugation with S). 
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